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B3.2.1  Introduction 

We  are  entering  an  era  when  condensed  matter  chemistry  arid  physics  can  be  predicted  from  theory  with 
increasing  realism  and  accuracy.  This  is  particularly  important  in  cases  where  experiments  lead  to  ambiguous 
conclusions,  for  regimes  in  which  there  still  exists  no  experimental  probe  and  for  predictions  of  the  properties 
of  modem  materials  in  order  to  select  the  most  promising  ones  for  synthesis  and  experimental  testing.  For 
example,  continuing  miniaturization  in  microelectronics  heightens  the  importance  of  understanding  of  quan¬ 
tum  effects,  which  computational  materials  theory  is  poised  to  provide,  based  to  some  degree  on  the  methods 
presented  here. 

Our  intention  is  to  give  a  brief  survey  of  advanced  theoretical  methods  used  to  determine  the  elec¬ 
tronic  and  geometric  structure  of  solids  and  surfaces.  The  electronic  structure  encompasses  the  energies  and 
wavefunctions  (and  other  properties  derived  from  them)  of  the  electronic  states  in  solids,  while  the  geometric 
structure  refers  to  the  equilibrium  atomic  positions.  Quantities  that  can  be  derived  from  the  electronic  struc¬ 
ture  calculations  include  the  electronic  (electron  energies,  charge  densities),  vibrational  (phonon  spectra), 
structural  (lattice  constants,  equilibrium  structures),  mechanical  (bulk  moduli,  elastic  constants)  and  opti¬ 
cal  (absorption,  transmission)  properties  of  crystals.  We  will  also  report  on  techniques  used  to  study  solid 
surfaces,  with  particular  examples  drawn  from  chemisorption  on  transition  metal  surfaces. 

In  his  chapter  on  the  fundamentals  of  quantum  mechanics  of  condensed  phases  (A1.3),  James  R  Che- 
likowsky  introduces  the  plane  wave  pseudopotential  method.  Here,  we  will  complement  his  chapter  by 
introducing  in  some  detail  tight-binding  methods  as  the  simplest  pedagogical  illustration  of  how  one  can 
construct  crystal  wavefunctions  from  atomic-like  orbitals.  These  techniques  are  very  fast  but  generally  not 
very  accurate.  After  reviewing  some  of  the  efforts  made  to  improve  upon  the  local  density  approximation 
(LDA,  explained  in  A  1.3),  we  will  discuss  general  features  of  the  technically  more  complex  all-electron  band 
structure  methods,  focusing  on  the  highly  accurate  but  not  very  fast  linear  augmented  plane  wave  (LAPW) 
technique  as  an  example.  We  will  introduce  the  idea  of  orbital-free  electronic  structure  methods  based  directly 
on  density  functional  theory  (DFT),  the  computational  effort  of  which  scales  linearly  with  size,  allowing  very 
large  systems  to  be  studied.  The  periodic  Hartree-Fock  (HF)  method  and  the  promising  quantum  Monte 
Carlo  (QMC)  techniques  will  be  briefly  sketched,  representing  many-particle  approaches  to  the  condensed 
phase  electronic  structure  problem. 

In  the  final  section,  we  will  survey  the  different  theoretical  approaches  for  the  treatment  of  adsorbed 
molecules  on  surfaces,  taking  the  chemisorption  on  transition  metal  surfaces,  a  particularly  difficult  to  treat 
yet  extremely  relevant  surface  problem  [1],  as  an  example.  While  solid  state  approaches  such  as  DFT  are  often 
used,  hybrid  methods, are  also  advantageous.  Of  particular  importance  in  this  area  is  the  idea  of  embedding, 
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aft  overview  of  Professor  Carter's  group's  work  using  pseudospectral  methods,  see: 

Martinez  T  J  and  Carter  E  A  1 995  Pseudospectral  methods  applied  to  the  electron  correlation  problem  Modern  Electronic  Struct 
Theory  vol  2,  ed  D  R  Yarkony  {Singapore:  World  Scientific)  pp  1 132-65  .  ure 

[75]  Hylleraas  E  A  and  Undheim  B  1 930  Z.  Phys.  65  759 

MacDonald  J  K  L  1933  Successive  approximations  by  the  Rayleigh-Ritz  variation  method  Phys.  Rev,  43  830-3 

[76]  Pople  J  A  1973  Theoretical  models  for  chemistry  Energy,  Structure,  and  Reactivitv  ed  D  W  Smith  and  W  B  McRae  (New  York- 

Wiley)  p  51-67 

[77]  Roos  B  0,  Taylor  P  R  and  Siegbahn  P  E  M  1 980  A  complete  active  space  SCF  method  (C ASSCF)  using  a  density  matrix  formulated 

super-CI  approach  Chem.  Phys.  48  157-73 

Roos  B  O  1987  The  complete  active  space  self-consistent  field  method  and  its  applications  in  electronic  structure  calculations 
Adv.  Chem.  Phys.  69  399-445 

[78]  -  Kelly  H  P  1963  Correlation  effects  in  atoms  Phyt.  Rev.  131  684-99 


[79]  Good  early  overviews  of  the  electron  propagator  (that  is  used  to  obtain  IP  and  EA  data)  and  of  the  polarization  propagator  are 

given  in:  ' 

Jorgensen  P  and  Simons  J  1981  Second  Quantization  Based  Methods  in  Quantum  Chemistry  (New  York:  Academic) 

The  very  early  efforts  on  these  methods  are  introduced  in: 

Linderberg  J  and  Ohm  Y  1973  Propagator  Methods  in  Quantum  Chemistry  (New  York:  Academic) 

More  recent  summaries  include: 

Cederbaum  L  S  and  Domcke  W  1977  Theoretical  aspects  of  ionization  potentials  and  photoelectron  spectroscopy  a  Green’s 
function  approach  Adv.  Chem.  Phys.  36  205-344 
Oddershede  J  1987  Propagator  methods  Adv.  Chem.  Phys.  69  201-39 

Ortiz  J  V  1997  The  electron  propagator  picture  of  moleculaj  electronic  structure  Computational  Chemistry;  Reviews  of  Current 
Trends  vol  2,  ed  J  Leszczynski  (Singapore:  World  Scientific)  pp  1-61 

[80]  The  introduction  of  EOMs  for  energy  differences  and  for  operators  that  connect  two  states  appears  first  in  the  nuclear  physics 

literature;  see  for  example: 

Rowe  D  J  1968  Equation-of-motion  method  and  the  extended  shell  model  Rev.  Mod.  Phys.  40  153-66 
I  applied  these  ideas  to  excitation  energies  in  atoms  and  molecules  in  1971;  see  equation  (2. 1  M2.6)  in: 

Simons  J  1971  Direct  calculation  of  first-  and  second-order  density  matrices.  The  higher  RPA  method  J.  Chem.  Phys.  55  1218-30 
In  1973,  the  EOM  method  was  then  extended  to  treat  IP  and  EA  cases: 

Simons  J  1973  Theory  of  electron  affinities  of  small  molecules  J.  Chem.  Phys.  58  4899-907 

In  a  subsequent  treatment  from  the  time-dependent  response  point  of  view,  connection  with  the  Greens  function  methods  was 
made: 

Simons  J  1972  Energy-shift  theory  of  low-lying  excited  electronic  states  of  molecules  J.  Chem.  Phys.  57  3787-92 
A  more  recent  overview  of  much  of  the  EOM,  Greens  function,  and  propagator  field  is  given  in: 

Oddershede  J  1987  Propagator  methods  Adv.  Chem.  Phys.  69  201-39 

[81]  Olsen  J  and  Jprgensen  P  1995  Time-dependent  response  theory  with  applications  to  self-consistent  field  and  multiconfigurational 

self-consistent  field  wave  functions  Modem  Electronic  Structure  Theory  vol  2,  ed  D  R  Yarkony  (Singapore:  World  Scientific) 
pp  857-990 

[82]  A  good  overview  of  the  recent  status  is  given  in: 

Bartlett  R  J  1995  Coupled-cluster  theory:  an  overview  of  recent  developments  Modern  Electronic  Structure  Theory  vol  2, 
ed  D  R  Yarkony  (Singapore:  World  Scientific)  pp  1047-131 

[83]  Helgaker  T,  Gauss  J,  J0rgensen  P  and  Olsen  J  1997  The  prediction  of  molecular  equilibrium  structures  by  the  standard  electronic 

wave  functions  J.  Chem.  Phys.  106  6430-40 
for  a  listing  and  for  further  details  on  this  study 

[84]  Two  review  papers  that  introduce  and  compare  the  myriad  of  semi-empirical  methods: 

Stewart  J  J  P  1991  Semiempirical  molecular  orbital  methods  Reviews  in  Computational  Chemistry  vol  1,  ed  K  B  Lipkowitz  and 
D  B  Boyd  (New  York:  VCH)  pp  45-8 1 

Zemer  M  C  1991  Semiempirical  molecular  orbital  methods  Reviews  in  Computational  Chemistry  vol  2,  ed  K  B  Lipkowitz  and 
D  B  Boyd  (New  York:  VCH)  313-65 

A  very  recent  overview,  including  efforts  to  interface  semi-empirical  electronic  structure  with  molecular  mechanics  treatments  of 
some  degrees  of  freedom  is  given  by: 

Thiel  W  1996  Perspectives  on  semiempirical  molecular  orbital  theory  New  Methods  in  Computational  Quantum  Mechanics  (Adv. 

Chem.  Phys.  XCIII)  ed  I  Prigogine  I  and  S  A  Rice  (New  York:  Wiley)  pp  703-57 
Earlier  texts  dealing  with  semi -empirical  methods  include: 

Pople  J  A  and  Beveridge  D  L  1910  Approximate  Molecular  Orbital  Theory  (New  York:  McGraw-Hill) 

Murrell  J  N,  Kettle  S  F  A  and  Tedder  J  M  1965  Valence  Theory  2nd  edn  (London:  Wiley) 
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where  a  small  cluster  of  surface  atoms  around  the  adsorbate  is  treated  with  more  care  than  the  surrounding 
region.  The  advantages  and  disadvantages  of  the  approaches  are  discussed.  : 

B3.2.2  Tight-binding  methods 

B 3. 2. 2.1  Tight  binding:  from  empirical  to  self-consistent 

The  wavefunction  in  a  solid  can  be  thought  to  originate  from  two  different  limiting  cases.  One  extreme  is  the 
nearly  free  electron  (NFE)  approach.  The  idea  here  is  that  the  valence  electrons  are  hardly  affected  by  the 
periodic  potential  of  the  atomic  cores.  Their  wavefunctions  can  then  be  assumed  to  be  easily  described  as 
linear  combinations  of  the  solutions  for  free  electrons:  the  plane  waves,  exp(i k  •  r).  The  NFE  approximation 
is  particularly  useful  for  so-called  NFE  metals,  such  as  the  alkali  metals.  At  the  other  extreme,  the  solid  can  be 
viewed  as  constructed  from  individual  atoms.  The  valence  wavefunctions  of  the  solid  are  then  approximated 
as  linear  combinations  of  the  wavefunctions  of  the  valence  electrons  of  the  atoms  (see  also  section  Al. 3.5.6). 
In  this  case,  the  electrons  are  considered  to  be  ‘tightly  bound’  to  the  atoms.  This  is  a  physically  reasonable 
view  of  covalently  bound  solids  and  molecules,  where  localized  chemical  bonds  are  the  norm  (bulk  silicon, 
organic  or  biomolecules  etc).  Methods  which  employ  this  view  of  the  electrons  in  the  solid  are  called  tight- 
binding  (TB)  methods.  The  wavefunctions  are  generally ^expanded  in  atomic  orbitals  (in  a  linear  combination 
of  atomic  orbitals  (LCAO)  formalism)  or  similarly  localized  functions. 

An  advantage  of  TB  is  that  generally  the  number  of  basis  functions  linearly  combined  to  give  the 
wavefunctions  is  rather  small.  The  solution  of  the  Schrodinger  equation  in  these  bases  is  then  fast  because 
the  matrices  representing  the  operators  are  small.  Also,  the  construction  of  the  Hamiltonian  matrix  elements 
is  fast,  since  generally  a  number  of,  sometimes  drastic,  approximations  are  made.  At  the  same  time,  however, 
the  small  basis  set  generally  limits  the  quality  of  the  TB  results,  since  the  variational  freedom  for  the  solution 
of  the  Schrodinger  equation  is  not  as  high  as  in  other  methods.  The  approximations  of  Hamiltonian  matrix 
elements  often  further  reduce  the  quality  of  the  results. 

Today,  the  term  TB  method  is  generally  understood  to  refer  to  a  technique  using  TB  basis  functions 
in  which  the  Hamiltonian  matrix  elements  are  adjusted  to  reproduce  results  from  experiments  and/or  from 
more  sophisticated  electronic  structure  methods  [2].  Depending  on  the  degree  of  dependence  on  external 
parameters,  the  methods  are  called  empirical  or  semi-empirical  TB.  A  number  of  approaches  are  used  for  the 
fitting  of  the  TB  parameters,  generally  a  tough  minimization  task  with  many  minima  (using  genetic  algorithms 
has  proved  quite  efficient  [3,4]).  It  has  been  noted  that  ‘great  care  is  needed  to  test  the  resulting  model  for 
reasonable  behavior  outside  the  range  of  the  fit’  [5, 6].  A  disadvantage  of  the  empirical  methods  is  that  it  is 
difficult  to  distinguish  to  what  extent  the  parametrization  or  the  method  itself  is  responsible  for  errors  in  the 
results. 

Frequent  approximations.made  in  TB  techniques  in  the  name  of  achieving  a  fast  method  are  the  use  of  a 
minimal  basis  set,  the  lack  of  a  self-consistent  charge  density,  the  fitting  of  matrix  elements  of  the  potential, 
the  assumption  of  an  orthogonal  overlap  matrix,  a  cut-off  radius  used  in  the  integration  to  determine  matrix 
elements,  and  the  neglect  of  matrix  elements  that  require  three-centre  integrals  and  crystal-field  terms.  We 
will  now  provide  more  details  on  these  approximations. 

Generally,  the  following  ansatz  for  the  wavefunction  is  made: 

ctl 

where  <pai(r)  =  {r\(pai)  represents  an  atomic  orbital  of  symmetry  a  (such  as  s,  p.v,  py,  p,)  at  atom  1. 

This  yields  the  generalized  eigenvalue  problem 

He'  =  €iSc\ 


(B3.2.1) 
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with  the  elements  of  the  Hamiltonian  matrix  Haipm  =  (<pai\H\  (ppm )  and  the  overlap  matrix  s  {<pai\<Ppm)- 

In  the  TB  approximation,  the  basis  functions  are  thought  to  be  sufficiently  localized  such  that  contributions  to 
the  Hamiltonian  matrix  usually  are  accounted  for  only  up  to  at  most  the  third  or  fourth  neighbour.  Frequently 
a  minimal  basis  set  is  used,  i.e.  a  single  orbital  <pai  is  used  per  atom  and  per  orbital  symmetry  to  expand  the 
wavefunction. 

In  orthogonal  TB  methods,  the  overlap  matrix  is  assumed  to  be  diagonal,  even  though  the  basis  functions 
of  adjacent  sites  ordinarily  are  not  orthogonal  [6],  Harrison  has  shown  that  this  approximation  can  be  compen¬ 
sated  for  by  adjustments  to  the  Hamiltonian  matrix  elements  (these  adjustments  are  arrived  at  automatically 
in  methods  depending  on  fitting,  for  example,  a  DFT  band  structure)  [7].  However,  this  approach  reduces 
the  transferability  of  the  TB  parameters  to  other  structures  [8],  Including  the  overlap  matrix  brings  with  it 
the  additional  cost  of  its  calculation  and  solving  the  generalized  eigenvalue  problem,  see  equation  (B3.2.1), 
rather  than  an  ordinary  eigenvalue  problem. 

One  can  construct  an  effective  potential,  written  here  in  the  DFT  language  (see,  for  example,  equa¬ 
tion  Al. 3.38  of  Al. 3)  as 

veff(r)  =  vext(r)  +  vH[/0(r)]  +  vxc[p(r)].  (B3.2.2) 

To  rationalize  the  ‘two-centre  approximation’,  the  effective  potential  is  written  as 

Ves(r)  =  2juefft/(|r  . 

/ 

where  ueff,/  is  centred  on  the  atom  /  and  vanishes  away  from  the  atom,  which  need  not  involve  any  approxi¬ 
mation. 

In  the  calculation  of  the  elements 


Halfim  —  {tPal 


T  +  ^  Veff,n 

n 


with  T  =  —  \  V2  the  kinetic  energy  operator,  several  types  of  potential  matrix  elements  can  be  distinguished 

[6]: 

(1)  Three-centre  terms,  i.e.  /  ^  m  ^  n.  These  are  frequently  neglected,  in  what  is  called  the  two-centre 
approximation,  based  on  the  assumed  strong  localization  of  the  orbitals  cpai(r ). 

(2)  Inter-atomic  two-centre  matrix  elements  {<pai\vtsj  +  vGfftm\<Ppm)-  These  matrix  elements  represent  the 
hopping  of  electrons  from  one  site  to  another.  They  can  be  described  [7]  as  linear  combinations  of 
so-called  Slater-Koster  elements  [9].  The  coefficients  depend  only  on  the  orientation  of  the  atoms  l 
and  m  in  the  crystal.  For  elementary  metals  described  with  s,  p,  and  d  basis  functions  there  are  ten 
independent  Slater-Koster  elements.  In  the  traditional  formulation,  the  orientation  is  neglected  and  the 
two-centre  elements  depend  only  on  the  distance  between  the  atoms  [6].  (In  several  models  [6, 10],  they 
have  been  made  dependent  on  the  environment  of  the  atoms  /  and  m.)  These  elements  are  generally  fitted 
to  reproduce  DFT  results  such  as  the  band  structure  or  the  values  of  DFT  matrix  elements  in  diatomics. 

(3)  Intra-atomic  matrix  elements,  or  on-site  terms,  with  l  =  m.  Traditionally,  the  potential  contributions 
from  other  atomic  sites,  ueff,n#/=m*  so-called  crystal-field  terms,  are  neglected  [10].  In  this  case,  then  the 
only  non-zero  on-site  terms  have  a  =  /?,  since  basis  functions  on  the  same  site  are  orthogonal  atomic 
orbitals.  There  are  methods  which  include  these  crystal-field  terms  [1 1, 12].  Physically,  these  diagonal 
elements  represent  the  energy  required  to  place  an  electron  in  a  specific  orbital.  In  some  implementations, 
they  are  set  to  the  orbital  energy  values  of  the  neutral  free  atom  [13],  guaranteeing  the  correct  limit  for 
isolated  atoms.  However,  this  approach  ignores  the  potential  contributions  to  the  diagonal  elements  due 
to  different  environments  in  a  molecule  or  crystal;  these  are  taken  into  account  in  other  variants  of  the 
method  [6, 10, 11]. 
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Most  TB  approaches  are  not  charge  self-consistent.  This  means  that  they  do  not  ensure  that  the  charge 
derived  from  the  wavefunctions  yields  the  effective  potential  vtff  assumed  in  their  calculation.  Some  methods 
have  been  developed  which  yield  charge  densities  consistent  with  the  electronic  potential  [14-16], 

The  localized  nature  of  the  atomic  basis  set  makes  it  possible  to  implement  a  linear-scaling  TB  algorithm, 
i.e.  a  TB  method  that  scales  linearly  with  the  number  of  electrons  simulated  [17].  (For  more  information  on 
linear  scaling  methods,  see  section  B3.2.3.3.) 

The  accuracy  of  most  TB  schemes  is  rather  low,  although  some  implementations  may  reach  the  accuracy 
of  more  advanced  self-consistent  LCAO  methods  (for  examples  of  the  latter  see  [18-20]).  However,  the 
advantages  of  TB  are  that  it  is  fast,  provides  at  least  approximate  electronic  properties  and  can  be  used  for 
quite  large  systems  (e.g.,  thousands  of  atoms),  unlike  some  of  the  more  accurate  condensed  matter  methods. 
TB  results  can  also  be  used  as  input  to  determine  other  properties  (e.g.,  photoemission  spectra)  for  which 
high  accuracy  is  not  essential. 

B3.2.2.2  Applications  of  tight-binding  methods 

TB  methods  have  been  widely  used  to  study  properties  of  simple  semiconductors  such  as  Si  [1 1]  and  GaN  [16]. 
In  the  latter  study,  the  effect  of  dislocations  on  the  electronic  structure  of  GaN  was  investigated  with  a  view 
toward  understanding  how  dislocations  affect  the  material’s  optical  properties.  The  large  supercell  of  224 
atoms  led  to  TB  as  the  method  of  choice.  This  particular  variant  of  TB  fits  TB  matrix  elements  to  DFT-LDA 
results  and  solves  self-consistently  for  atomic  charges.  It  has  also  been  used  to  predict  reaction  energetics  of 
organic  molecules,  the  structure  of  large  biomolecules  and  the  surface  geometry  and  band  structure  of  III-V 
semiconductors  [15].  The  TB  method  is  expected  to  provide  qualitatively  reasonable  results  for  systems 
where  localized  atomic  charges  make  sense  and  hence  is  not  expected  to  perform  as  well  for  metallic  systems. 
Despite  potential  problems  of  TB  for  metals,  the  TB  approach  has  also  been  used  to  study  the  phonon  spectrum 
of  the  transition  metal  molybdenum  [6],  the  elastic  constants,  vacancies  and  surfaces  of  monatomic  transition 
and  noble  metals  and  the  Hall  coefficient  of  complex  perovskite  crystals  [10].  As  an  example  of  data  available 
from  a  TB  calculation,  a  TB  variant  of  extended  Hiickel  theory  [2 1 , 22]  was  used  to  describe  the  initial  states  in 
photoemission  from  GaN  [23].  The  parameters  were  fitted  to  the  bulk  band  structure  E„ (fc)  (for  a  definition, 
see  section  Al.3.6).  As  displayed  in  figure  B3.2.1,  good  agreement  is  found  for  the  occupied  states  (negative 
energies),  while  larger  differences  for  the  conduction  bands  (positive  energies)  reveal  a  typical  problem  of 
the  TB  methods:  they  are  far  less  capable  of  describing  the  delocalized  conduction  band  states  (the  same  is 
true  for  delocalized  valence  states  in  a  metal,  as  mentioned  above).  In  figure  B3.2.2,'we  show  a  series  of 
calculated  photoemission  spectra  compared  to  experimental  results  [23].  The  dispersion  of  the  main  peaks  as 
a  function  of  emission  angle  and  photon  energy  agrees  reasonably  well  in  theory  and  experiment. 

B3.2.3  First-principles  electronic  structure  methods 

In  this  section,  we  briefly  review  the  basic  elements  of  DFT  and  the  LDA.  We  then  focus  on  improvements 
suggested  to  remedy  some  of  the  shortcomings  of  the  LDA  (see  section  B  3. 2. 3. 1).  A  wide  variety  of  techniques 
based  on  DFT  have  been  developed  to  calculate  the  electron  density.  Many  approaches  do  not  calculate  the 
density  directly  but  rather  solve  for  either  a  set  of  single-electron  orbitals,  or  the  Green’s  function,  from  which 
the  density  is  derived. 

In  section  B3.2.3.2,  we  introduce  a  number  of  techniques  commonly  referred  to  as  ab  mirier  all-electron 
electronic  structure  methods.  Ab  initio  methods,  in  particular,  aim  at  calculating  the  energies  of  electrons 
and  their  wavefunctions  as  accurately  as  possible,  introducing  as  few  adjustable  parameters  as  possible. 
(Empirical  or  semi-empirical  methods  include  the  empirical  pseudopotential  approach  (see  section  Al.3.5.5) 
and  many  TB  techniques  (see  section  B3.2.2).)  Within  the  ab  initio  band  structure  approach,  two  communities 
exist  that  differ  in  their  treatment  of  the  singular  nature  of  realistic,  Coulomb-like  crystal  potentials.  In  the 
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Figure  B3.2.1.  The  band  structure  of  hexagonal  GaN,  calculated  using  EHT-TB  parameters  determined  by  a  genetic 
algorithm  [23].  The  target  energies  are  indicated  by  crosses.  The  target  band  structure  has  been  calculated  with  an 
ab  initio  pseudopotential  method  using  a  quasiparticle  approach  to  include  many-particle  corrections  [14], 


pseudopotential  approach  discussed  by  Chelikowsky  in  chapter  A1.3,  the  Coulomb  singularity  (  Z/r)  o 
the  crystal  potential  is  replaced  by  a  smoother  function,  whereas  in  the  so-called  ‘all-electron  approach,  the 
Coulomb  singularity  is  retained.  The  pseudopotential  transformation  limits  the  range  of  electron  energies 
which  can  be  accessed.  However,  since  the  pseudo-wavefunction  is  much  smoother  than  the  all-electron 
wavefunction  (which  has  large  oscillations  near  the  nucleus),  the  pseudopotential  allows  the  use  of  a  plane 
wave  basis  set,  which  is  comparatively  easy  to  handle.  In  principle,  the  all-electron  methods  have  no  hmi  a  ion 
on  the  energy  range  of  calculations.  This  is  achieved  by  a  sophisticated  representation  of  the  wavefunction 
The  so-called  orbital-free  DFT  technique,  which  aims  to  directly  calculate  the  electron  density  for  which 
the  total  energy  is  minimal,  is  presented  as  an  example  of  methods  whose  computational  effort  scales  linear  y 
with  system  size  (see  section  B3.2.3.3).  In  section  B3.2.3.4,  we  discuss  the  periodic  HF  method,  an  alternative 
approach  to  DFT  that  offers  a  well  defined  starting  point  for  many-particle  corrections.  Finally,  the  two  most 

frequently  used  QMC  techniques  are  described  in  section  B3.2.3.5. 

B3.2.3.1  The  local  density  approximation  and  beyond 

In  DFT,  the  electronic  density  rather  than  the  wavefunction  is  the  basic  variable.  Hohenberg  and  Kohn  showed 
f24]  that  all  the  observable  ground-state  properties  of  a  system  of  interacting  electrons  moving  in  an  externa 
potential  vM(r)  are  uniquely  dependent  on  the  charge  density  p(r)  that  minimizes  the  system  s  total  energy. 
However  there  is  no  known  formula  to  calculate  from  the  density  the  total  energy  of  many  electrons  moving 
in  a  general  potential.  Hohenberg  and  Kohn  proved  that  there  exists  a  universal  functional  of  the  density, 
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Figure  B3.2.2.  A  series  of  photoemission  spectra.  The  angles  give  the  polar  angle  of  electron  emission  at  the  stated  photon 
energy  scanning  the  surface  Brillouin  zone  from  f  to  M.  Left:  A  calculation  using  the  tight-binding  parametrization 
(given  the  band  structure  in  figure  B3.2. 1 )  for  the  initial  states  [23].  Right:  Experimental  spectra  by  Dhesi  et  al  [195],  The 
difference  in  binding  energies  is  due  to  the  experimental  difficulty  in  determining  the  Fermi  energy  [23].  (Experimental 
figure  by  Professor  K  E  Smith.) 


called  G[p],  such  that  the  expression 


>  E[p]  =  J  uext(r)p(r)  dV  +  i  j  ^^pd3 r  dV  +  G[p]  (B3.2.3) 

has  as  its  minimum  value  the  correct  ground-state  energy  associated  with  uext(r)-  Here,  the  first  term  on  the 
right-hand  side  represents  the  energy  due  to  an  external  potential,  including  the  electron-nuclear  potential, 
while  the  second  term  is  the  classical  Coulomb  energy  of  the  electronic  system.  The  functional  G[p]  is  valid 
for  any  number  of  electrons  and  any  external  potential,  but  it  is  unknown  and  further  steps  are  necessary  to 
approximate  it. 

Kohn  and  Sham  [25]  decompose  G[p]  into  the  kinetic  energy  of  an  analogous  set  of  non-interacting 
electrons  with  the  same  density  p(r)  as  the  interacting  system, 


i  ' 


(where  ^(r)  (rl^-)  is  the  wavefunction  of  electron  /),  and  the  exchange  and  correlation  energy  of  an 

interacting  system  with  density  p(r),  Exc[p].  The  functional  £xc[/>]  is  not  known  exactly.  Physically,  it 
represents  all  the  energy  corrections  beyond  the  Hartree  term  to  the  independent-particle  model,  i.e.  the  non- 
classical  many-body  effects  of  exchange  and  correlation  (xc)  and  the  difference  between  the  kinetic  energy 
of  the  interacting  electron  system  T[p]  and  the  analogous  non-interacting  system  Ts[p]. 
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In  the  LDA,  the  exchange  and  correlation  energy  is  approximated  using  the  exchange  and  correlation 
energy  of  the  homogeneous  electron  gas  at  the  same  density  (see  section  Al.3.3.3).  The  crystal  density  is 
obtained  by  solving  the  single-particle  Kohn-Sham  equation 

■  V"  +  ueff(r)  j  fj{r)  =  Eifi(r),  (B3.2.4) 

for  a  self-consistent  potential  veff,  i.e.  a  potential  which  is  produced  by  the  density  p.  In  bulk  crystal  calcu¬ 
lations,  the  index  i  runs  over  both  the  Bloch  vector  fc  (see  section  A  1.3.4)  and  the  band  index  n  (in  a  simple 
crystal,  this  band  could  be  derived,  for  example,  entirely  from  s  states).  The  solutions  to  equation  (B3.2.4) 
are  often  called  Kohn-Sham  orbitals.  The  crystal  density  is  then 

p(r)  = 

I  ■ 

The  eigenenergy  £,  can  be  defined  as  the  derivative  of  the  total  energy  of  the  many-electron  system  with 
respect  to  the  occupation  number  of  a  specific  orbital  [26].  In  HF  theory  (where  equation  (B3.2.4)  applies  and 
the  ueff  contains  a  non-local  exchange  operator,  see  section  Al.3.1.2  and  chapter  B3.1),  Koopmans  theorem 
states  that  the  single-particle  eigenvalue  is  the  negative  of  the  ionization  energy  (neglecting  the  relaxation  of 
the  electronic  system).  In  contrast,  the  identification  of  the  highest  occupied  Kohn-Sham  eigenvalue  with 
the  negative  of  the  ionization  energy  is  a  controversial  subject  [27],  While  there  is  no  rigorous  connection 
between  eigenvalue  differences  and  excitation  energies  in  either  HF  or  DF  theory,  comparisons  of  these  values 
are  common  practice  (see  below  for  more  appropriate  methods).  Relative  differences  among  occupied  single¬ 
particle  energies  often  agree  well  with  the  experiment.  Even  though  DFT  only  provides  a  solution  for  the 
ground  state  of  the  electronic  system,  the  energy  differences  in  the  lower  conduction  bands,  i.e.  low-energy 
excited  states,  often  are  represented  surprisingly  well,  too.  However,  in  LDA  calculations  of  semiconductors 
and  insulators,  almost  always  the  size  of  the  gap  between  the  valence  band  maximum  and  the  conduction 
band  minimum  is  underestimated,  since  many-particle  effects  are  incorrectly  represented  by  the  parametrized 
exchange-correlation  energy  (see,  for  example,  [28]).  One  ad  hoc  remedy,  which  works  well  for  many  systems 
and  which  is  employed  in  the  examples  presented  here,  is  to  use  what  is  amusingly  referred  to  as  a  scissor 
operator ,  i.e.  a  rigid  shift,  to  correct  the  gap  size  [29, 30].  Typically  the  shift  is  determined  by  knowing,  for 
example,  the  DFT  error  in  predicting  the  measured  optical  band  gap.  The  entire  conduction  band  is  shifted 
rigidly  upward  by  the  amount  to  match  the  experimental  band  gap. 

More  advanced  techniques  take  into  account  quasiparticle  corrections  to  the  DFT-LDA  eigenvalues. 
Quasiparticles  are  a  way  of  conceptualizing  the  elementary  excitations  in  electronic  systems.  They  can  be 
determined  in  band  structure  calculations  that  properly  include  the  effects  of  exchange  and  correlation.  In 
the  LDA,  these  effects  are  modelled  by  the  exchange-correlation  potential  i4DA.  In  order  to  more  accurately 
account  for  the  interaction  between  a  particle  and  the  rest  of  the  system,  the  notion  of  a  local  potential  has 
to  be  generalized  and  a  non-local,  complex  and  energy-dependent  exchange-correlation  potential  has  to  be 
introduced,  referred  to  as  the  self-energy  operator  E(t\  r';  E).  The  self-energy  can  be  expanded  in  terms  of 
the  screened  Coulomb  potential  W ,  where  W  =  e~lv  is  the  Coulomb  interaction  v  screened  by  the  inverse 
dielectric  function  e-1 .  In  a  lowest  order  expansion  in  W,  the  self-energy  can  be  approximated  as  E  =  GW, 
giving  the  GW  approximation  [31].  Here  G  is  the  one-electron  Green’s  function  describing  the  propagation 
of  an  additional  electron  injected  into  a  system  of  other  electrons  (it  can  also  describe  the  extraction  of  an 

electron).  . 

To  be  a  bit  more  explicit  (following  [32, 33]),  the  quasiparticle  energies  and  wavefunctions  are  given  by 

( T  +  vext  +  uh  +  f  dr'E(r,  r;  E„k)irnk(r')  =  Enk^nkir), 


where  T  is  the  kinetic  energy  operator,  uext  is  the  external  potential  due  to  the  ions,  and  uH  is  the  Hartree 
Coulomb  interaction.  Since  the  self-energy  operator  in  general  is  non-Hermitian,  the  quasiparticle  energies 
Enk  are  complex  in  general,  and  the  imaginary  part  gives  the  lifetime  of  the  quasiparticle.  To  first  order  in 
W ,  the  self-energy  is  then  given  by  •  ; 

E(r,  r'\  E)  =  J-  f  da)e~'SulG(r,  r';  E  -  <o)W(r,  r'\  to)  ‘ 

2  it  J  _ 


where  8  is  a  positive  infinitesimal  and  cd  corresponds  to  an  excitation  frequency.  The  inputs  are  the  full 
interacting  Green’s  function, 

G(r,/;£)  =  v|a(r) 

E~  E<ik  -  iS,,k 

where  is  an  infinitesimal  and  the  dynamically  screened  Coulomb  interaction. 


W(r,  r\  co)  =  Q  /  d rf,€  1  (r,  rft;  co)v(rfi 


where  is  the  inverse  dielectric  matrix,  v(r)  =  l/|r|  and  Q  is  the  volume  of  the  system.  Usually  the 
calculations  start  with  the  construction  of  the  Green’s  function  and  the  screened  Coulomb  potential  from 
self-consistent  LDA  results.  The  self-energy  2  then  has  to  be  obtained  together  with  G  in  a  self-consistent 
procedure.  However,  due  to  the  severe  computational  cost  of  this  procedure,  it  is  usually  not  carried  out  (see, 
for  example,  [34]).  Instead,  it  is  common  practice  to  construct  the  self-energy  operator  non-self-consistently 
using  the  self-consistent  LDA  results  to  determine  quasiparticle  corrections  to  the  LDA  energies,  resulting 
in  the  quasiparticle  band  structure.  The  GW  approximation  has  been  applied  to  a  wide  range  of  metals, 
semiconductors  and  insulators,  where  it  has  been  found  to  lead  to  striking  improvements  in  the  agreement  of 
optical  excitation  spectra  with  the  experiment  (see,  for  example  [32, 35-37]).  Recent  studies  also  found  that 
the  GW  charge  density  is  close  to  the  experiment  for  diamond  structure  semiconductors  [38],  and  lifetimes 
of  low-energy  electrons  in  metals  have  been  calculated  [39]. 

Another  disadvantage  of  the  LDA  is  that  the  Hartree  Coulomb  potential  includes  interactions  of  each 
electron  with  itself,  and  the  spurious  term  is  not  cancelled  exactly  by  the  LDA  self-exchange  energy,  in 
contrast  to  the  HF  method  (see  A  1.3),  where  the  self-interaction  is  cancelled  exactly.  Perdew  and  Zunger 
proposed  methods  to  evaluate  the  self-interaction  correction  (SIC)  for  any  energy  density  functional  [40]. 
However,  full  SIC  calculations  for  solids  are  extremely  complicated  (see,  for  -example*  [41-43]).  As  an 
alternative  to  the  very  expensive  GW  calculations,  Pollmann  etal  have  developed  a  pseudopotential  built  with 
self- interaction  and  relaxation  corrections  (SIRC)  [44].  The  pseudopotential  is  derived  from  an  all-electron 
SIC-LDA  atomic  potential.  The  relaxation  correction  takes  into  account  the  relaxation  of  the  electronic 
system  upon  the  excitation  of  an  electron  [44].  The  authors  speculate  that  \  .  .the  ability  of  the  SIRC  potential 
to  produce  considerably  better  band  structures  than  DFT-LDA  may  reflect  an  extra  nonlocality  in  the  SIRC 
pseudopotential,  related  to  the  nonlocality  Or  orbital  dependence  in  the  SIC  all-electron  potential.  In  addition, 
it  may  mimic  some  of  the  energy  and  the  non-local  space  dependence  of  the  self-energy  operator  occurring 
in  the  GW  approximation  of  the  electronic  many  body  problem’ [45]. 

The  LDA  also  fails  for  strongly  correlated  electronic  systems.  Examples  of  such  systems  are  the  late 
3d  transition-metal  mono-oxides  MnO,  FeO,  CoO,  and  NiO.  Within  the  local  spin  density  approximation 
(LSDA),  the  energy  gaps  calculated  for  MnO  and  NiO  are  too  small  [46]  and,  even  worse,  FeO  and  CoO 
are  predicted  to  be  metallic,  whereas  experimentally  they  have  been  found  to  be  large-gap  insulators.  While 
the  GW  approximation  yields  an  energy  gap  of  NiO  in  reasonable  agreement  with  experiment  [47],  the 
computational  cost  of  this  procedure  is  very  high.  The  SIC-LDA  method  reproduces  quite  well  the  strong 
localization  of  the  d  electrons  in  transition  metal  compounds,  but  the  orbital  energies  obtained  by  SIC  are 
usually  in  strong  disagreement  with  experimental  results  (for  transition  metal  oxides,  for  example,  occupied 
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d  bands  are  approximately  \  Hartree  below  the  oxygen  valence  band — a  separation  not  seen  in  spectroscopic 
data:  see,  for  example,  the  experimental  results  in  [48])  [49).  An  alternative  solution  to  this  problem  is  offered 
by  the  LDA +U  method  [49, 50],  where  LDA  encompasses  the  LSDA.  In  the  LDA+t/  technique,  the  electrons 
are  divided  into  two  subsystems  which  are  treated  separately:  the  strongly  localized  (d  or  f)  electrons  and 
the  delocalized  s  and  p  electrons.  The  latter  are  treated  by  standard  LDA.  The  on-site  interactions  among  the 
strongly  localized  electrons  on  each  atom,  however,  are  taken  into  account  by  a  term  nini'  where 

n.  are  the  occupation  numbers  of  the  strongly  localized  orbitals  and  U  is  the  Coulomb  interaction  parameter 
(for  details  on  the  first-principles  calculation  of  U,  see  [51]).  At  least  for  localized  d  or  f  states,  the  LDA+t/ 
technique  may  be  viewed  as  an  approximation  to  the  GW  approximation  [49].  Band  gaps,  valence  band 
widths  and  magnetic  moments  have  been  calculated  with  LDA+t/  that  agree  with  experiment  for  a  variety  of 
transition  metal  compounds  [49, 52],  among  other  applications. 

J 33  2.3.2  All-electron  DFT  methods 

(a)  Introduction 

When  the  highest  accuracy  is  sought  for  the  electronic  and  geometric  properties  of  crystals,  all  the  electrons 
of  the  atoms  in  the  crystal  and  the  full  Coulomb  singularity  of  the  nuclear  potential  must  be  accounted  for. 
All-electron  approaches,  which  do  just  that,  generally  cannot  compete  with  pseudopotential  techniques  in 
speed  and  simplicity  of  algorithm.  However,  the  latter  suffer  from  severe  drawbacks  when  it  comes  to  the 
construction  of  the  very  pseudopotentials  these  methods  depend  upon:  even  for  so-called  ab  initio  potentials, 
the  pseudopotentials  are  far  from  uniquely  determined.  Additionally/problems  with  transferability  and  the 
construction  of  potentials  for  such  elements  as  the  transition  metals  remain.  All-electron  techniques  can  deal 
with  any  element  and  there  are  no  worries  about  transferability  of  the  potential.  However,  the  accuracy  comes 
at  a  price:  due  to  the  Coulomb  singularity  of  the  potential  at  the  nuclear  positions,  the  wavefunctions  are 
highly  oscillatory  close  to  the  nucleus.  For  those  all-electron  methods  that  use  wavefunctions  to  represent 
the  electrons  (a  Green’s  function  method,  for  example,  does  not),  this  means  that  a  simple  plane  wave  basis 
set  cannot  be  used  for  the  expansion  of  the  wavefunctions.  To  reach  convergence  of  a  plane  wave  exp(i/c  ■  r) 
expansion  would  require  a  prohibitive  number  of  basis  functions.  Thus,  specialized  basis  sets  have  been 
invented  for  all-electron  calculations. 

We  now  discuss  the  most  important  theoretical  methods  developed  thus  far:  the  augmented  plane  wave 
(APW)  and  the  Korringa-Kohn-Rostoker  (KKR)  methods,  as  well  as  the  linear  methods  (linear  APW  (LAPW), 
the  linear  muffin-tin  orbital  [LMTO]  and  the  projector-augmented  wave  [PAW])  methods. 

In  the  early  all-electron  techniques,  the  crystal  was  separated  into  spheres  around  the  atoms,  so-called 
‘muffin- tin’  spheres,  and  the  interstitial  region  in  between.  Inside  the  spheres,  the  potential  was  approximated 
as  spherically  symmetric,  while  in  the  interstitial  region  it  was  assumed  to  be  constant.  This  shape  approx¬ 
imation  of  the  potential  is  reasonable  for  close-packed  crystals  such  as  hexagonally  close-packed  metals, 
where  the  spheres  cover  a  large  fraction  of  the  crystal  volume.  However,  in  less  densely  arranged  crystals, 
such  as  diamond  structure  semiconductors  (see  figure  Al.3.4),  the  muffin-tin  approximation  leads  to  large 
errors.  In  the  diamond  and  the  related  zincblende  structures,  only  34%  of  the  volume  is  covered  by  touching 
muffin-tin  spheres  (figure  B3.2.3).  For  all  of  the  all-electron  methods,  versions  have  been  developed  that  are 
not  restricted  to  shape  approximations  of  the  potential.  These  techniques  are  referred  to  as  general,  or  full, 
potential  methods.  . 

(b)  The  augmented  plane  wave  method 

The  APW  technique  was  proposed  by  Slater  in  1937  [53,54].  It  remains  the  most  accurate  of  the  band  structure 
methods  for  the  muffin-tin  approximation  of  the  potential.  The  wavefunction  is  expanded  in  basis  functions 
(ft  (fc  +  Qh  Et  r),  the  APWs,  each  of  which  is  identical  to  the  plane  wave  exp(i(fc  +  G-t)  ■  r)  in  the  interstitial 
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equation  for  a  given  gy-  inr  /nnn  cincmlar')  at  the  origin.  With  the  basis  functions  cpdk  +  G,-,  E,  r), 
sphere  boundary  and  must  be  regular  (no  -  0  °  .  Since  the  Hamiltonian  matrix 

secular  equation  is  solved  by  finding  the  roots  of  the  determinant  of  the  H£E)  —  £S(£)  matnx.  P 

cannot  be  treated  by  the  eigenvalue  routines  of  linear  algebra.) 
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so-called ‘asymptote  problem  ).  . 
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1967,  the  LAPW  provided  a  way  to  avoid  the  above-mentioned  drawbacks  of  the  APW  technique, 
discuss. 
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Figure  B3.2.4.  A  schematic  illustration  of  an  energy-independent  augmented  plane  wave  basis  function  used  in  the 
LAPW  method.  The  black  sine  function  represents  the  plane  wave,  the  localized  oscillations  represent  the  augmentation 
of  the  function  inside  the  atomic  spheres  used  for  the  solution  of  the  Schrodinger  equation.  The  nuclei  are  represented  by 
filled  black  circles.  In  the  lower  part  of  the  picture,  the  crystal  potential  is  sketched. 


( c )  The  linear  augmented  plane  wave  method 

The  main  disadvantage  of  the  APW  technique  is  that  it  leads  to  a  nonlinear  secular  problem  because  the  basis 
functions  depend  on  the  energy.  A  number  of  attempts  have  been  made  to  construct  linear  versions  of  the  APW 
approach  by  introducing  energy-independent  basis  functions  in  different  ways.  In  1970,  Koelling  invented 
the  alternative  APW  [56]  and  Bross  the  modified  APW  [57].  In  1975,  Andersen  constructed  the  LAPW  [58] 
formalism,  which  today  is  the  most  popular  APW-like  band  structure  method.  Further  extensions  of  the  linear 
methods  appeared  in  the  early  1990s:  Singh  developed  the  LAPW  plus  localized  orbitals  (LAPW+LO  [59]) 
in  1991  and  Krasovskii  the  extended  LAPW  (ELAPW  [60])  in  1994.  Recently  the  APW+LO  technique  has 
been  implemented  by  Sjostedt  and  Nordstrom  [61]  according  to  an  idea  by  Singh.  While  the  LAPW  technique 
is  generally  used  in  combination  with  DFT  approaches,  it  has  also  been  applied  based  on  the  LDA+f/  [62] 
and  HF  theories  [63]. 

The  LAPW  method,  as  suggested  in  1975  [58,64],  avoids  the  problem  of  the  energy  dependence  of 
the  Hamiltonian  matrix  by  introducing  energy-independent  APW  basis  functions.  Here,  too,  the  APWs  are 
derived  from  plane  waves  by  augmentation:  Bessel  functions  fi{\ k  +  Gi\r)  in  the  Rayleigh  decomposition 
inside  the  muffin-tin  sphere  are  replaced  by  functions  uu(r)  derived  from  the  spherical  potential,  which 
are.  independent  of  the  energy  of  the  state  that  is  sought  and  that  match  the  Bessel  functions  at  the  sphere 
radius  in  value  and  in  slope  (see  figure  B3.2.4).  The  plane  wave  part  of  the  basis  remains  the  same  but  the 
energy-independent  APWs  allow  the  energies  and  the  wavefunctions  to  be  determined  by  solving  a  standard 
generalized  eigenvalue  problem. 

In  linearizing  the  APW  problem  as  it  is  done  in  the  LAPW  method,  the  variational  freedom  of  the  APW 
basis  set  is  reduced.  The  reason  is  that  the  wavefunction  inside  the  spheres  is  rigidly  coupled  to  its  plane  wave 
expansion  in  the  interstitial  region  [65].  This  means  that  the  method  cannot  yield  an  accurate  wavefunction 
even  if  the  eigenvalue  is  within  a  few  eV  of  the  chosen  energy  parameters  [66].  Flexibility  is  defined  in  this 
context  as  the  possibility  to  change  the  wavefunction  inside  the  spheres  independently  from  the  wavefunction 
in  the  interstitial  region.  Flexibility  can  be  achieved  in  the  linear  band  structure  methods  by  adding  basis 
functions  localized  inside  the  spheres  whose  value  and  slope  vanish  at  the  sphere  boundary  [54, 67, 68].  A 
‘flexible’  basis  set  extending  the  LAPW  with  localized  functions  is  preferable  to  the  one  used  in  the  pure 
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LAPW  technique.  Flexible  linear  methods  are  the  MAPW,  the  LAPW+LO  and  the  ELAPW,  the  latter  of 
which  provides  a  necessary  degree  of  flexibility  with  a  minimal  number  of  basis  functions  [65]. 

The  additional  functions  increase  the  matrix  dimension  slightly  and  thus  the  computational  effort.  How¬ 
ever,  the  increased  flexibility  of  the  basis  set  makes  possible  a  number  of  extensions  of  the  LAPW  method. 
One  is  a  k-p  formulation  of  the  ELAPW  method  [68],  which  would  lead  to  large  errors  in  the  regular  LAPW 
due  to  its  lesser  flexibility.  The  augmented  Fourier  components  (AFC)  technique  [69]  for  treating  a  general 
potential  is  based  on  this.  The  AFC  method  is  an  alternative  to  the  full-potential  LAPW  (FLAPW)  method 
[70, 71],  (Recently  progress  has  been  made  in  increasing  the  computational  efficiency  of  the  FLAPW  method 
[72].)  The  AFC  method  does  not  have  the  same  demanding  convergence  criteria  as  the  FLAPW  method  but 
yields  physically  equivalent  results  [69]. 

The  general  potential  LAPW  techniques  are  generally  acknowledged  to  represent  the  state  of  the  art  with 
respect  to  accuracy  in  condensed  matter  electronic-structure  calculations  (see,  for  example,  [62, 73]).  These 
methods  can  provide  the  best  possible  answer  within  DFT  with  regard  to  energies  and  wavefunctions. 

(d)  The  Korringa-Kohn-Rostoker  technique 

The  KKR  method  uses  multiple-scattering  theory  to  solve  the  Kohn-Sham  equations  [74, 75].  Rather  than 
calculate  the  wavefunction,  modem  incarnations  calculate  the  Green’s  function  G.  The  Green’s  function  is 
the  solution  to  the  equation  schematically  given  by  (H  -  E)G{E )  =  -<$,  where  H  is  the  Hamiltonian,  E  the 
single-electron  energy  and  8  the  delta  function  8(r-  r').  The  properties  of  the  system,  such  as  the  electron 
density,  the  density  of  states  and  the  total  energy  can  be  derived  from  the  Green’s  function  [73],  The  crystal  is 
represented  as  a  sum  of  non-overlapping  potentials;  in  the  modem  version,  there  are  no  shape  approximations, 
i.e.  the  potentials  are  space-filling  [76].  Within  the  multiple-scattering  formalism,  the  wavefunction  is  built 
up  by  taking  into  account  the  scattering  and  rescattering  of  a  free-electron  wavefunction  by  scatterers.  The 
scatterers  are  (generally)  the  atoms  of  the  crystal  and  the  single-scattering  properties  (the  properties  of  the 
isolated  scatterer)  are  derived  from  the  effective,  singular  potentials  of  the  atoms  (given  in  equation  (B3.2.2)). 
The  Green’s  matrix  is  then  constructed  from  the  knowledge  of  the  scattering  properties  of  the  single  scatterers 
and  the  analytically  known  Green’s  function  of  the  free  electron.  The  full-potential  KKR  method  has  been 
shown  to  have  the  same  level  of  accuracy  as  the  full-potential  LAPW  method  [73],  The  Green’s  function 
formulation  offers  the  advantage  of  easy  inclusion  of  defects  in  the  bulk  or  clean  surfaces.  Such  calculations 
start  with  the  Green’s  function  of  the  periodic  crystal  and  include  the  perturbation  through  a  Dyson  equation 
[77],  Yussouff  states  that  the  difference  in  speeds  between  the  linear  methods  and  his  ‘fast’  KKR  technique 
is  at  most  a  factor  of  ten,  in  favour  of  the  former  [78].  While  the  KKR  technique  has  an  accuracy  comparable 
to.the  APW  method,  it  has  the  disadvantage  of  not  being  a  linear  approach,  limiting  speed  and  simplicity. 

(e)  The  linear  muffin-tin  orbital  method 

The  LMTO  method  [58, 79]  can  be  considered  to  be  the  linear  version  of  the  KKR  technique.  According 
to  official  LMTO  historians,  the  method  has  now  reached  its  ‘third  generation’  [79]:  the  first  starting  with 
Andersen  in  1975  [58],  the  second  commonly  known  as  TB-LMTO.  In  the  LMTO. approach,  the  wavefunction 
is  expanded  in  a  basis  of  so-called  muffin-tin  orbitals.  These  orbitals  are  adapted  to  the  potential  by  constructing 
them  from  solutions  of  the  radial  Schrodinger  equation  so  as  to  form  a  minimal  basis  set.  Interstitial  properties 
are  represented  by  Hankel  functions,  which  means  that,  in  contrast  to  the  LAPW  technique,  the  orbitals  are 
localized  in  real  space.  The  small  basis  set  makes  the  method  fast  computationally,  yet  at  the  same  time  it  . 
restricts  the  accuracy.  The  localization  of  the  basis  functions  diminishes  the  quality  of  the  description  of  the 
wavefunction  in  the  interstitial  region. 

In  the  commonly  used  atomic  sphere  approximation  (ASA)  [79],  the  density  and  the  potential  of  the 
crystal  are  approximated  as  spherically  symmetric  within  overlapping  muffin-tin  spheres.  Additionally,  all 
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integrals,  such  as  for  the  Coulomb  potential,  are  performed  only  over  the  spheres.  The  limits  on  the  accuracy 
of  the  method  imposed  by  the  ASA  can  be  overcome  with  the  full-potential  version  of  the  LMTO  (FP-LMTO) 
which  gives  highly  accurate  total  energies  [79, 80].  It  was  found  that  the  FP-LMTO  is  ‘at  least  as  accurate  as, 
and  much  faster  than,’  pseudopotential  plane  wave  calculations  in  the  determination  of  structural  and  dynamic 
properties  of  silicon  [80].  The  FP-LMTO  is  considerably  slower  than  LMTO- AS  A,  however,  and  it  has  been 
found  that  ASA  calculations  can  yield  accurate  results  if  the  full  expansion,  rather  than  only  the  spherical 
part,  of  the  charge  is  used  in  what  is  called  a  full-charge  (rather  than  a  full-potential)  method  and  the  integrals 
are  performed  exactly  [73, 79]. 

The  LMTO  method  is  the  fastest  among  the  all-electron  methods  mentioned  here  due  to  the  small  basis 
size.  The  accuracy  of  the  general  potential  technique  can  be  high,  but  LAPW  results  remain  the  ‘gold  standard’ . 

(f)  The  projector  augmented  wave  technique 

The  projector  augmented-wave  (PAW)  DFT  method  was  invented  by  Blochl  to  generalize  both  the  pseudopo¬ 
tential  and  the  LAPW  DFT  techniques  [81].  PAW,  however,  provides  all-electron  one-particle  wavefunctions 
not  accessible  with  the  pseudopotential  approach.  The  central  idea  of  the  PAW  is  to  express  the  all-electron 
quantities  in  terms  of  a  pseudo-wavefunction  (easily  expanded  in  plane  waves)  term  that  describes  interstitial 
contributions  well,  and  one-centre  corrections  expanded  in  terms  of  atom-centred  functions,  that  allow  for 
the  recovery  of  the  all-electron  quantities.  The  LAPW  method  is  a  special  case  of  the  PAW  method  and 
the  pseudopotendal  formalism  is  obtained  by  an  approximation.  Comparisons  of  the  PAW  method  to  other 
all-electron  methods  show  an  accuracy  similar  to  the  FLAPW  results  and  an  efficiency  comparable  to  plane 
wave  pseudopotential  calculations  [82, 83].  PAW  is  also  formulated  to  carry  out  DFT  dynamics,  where  the 
forces  on  nuclei  and  wavefunctions  are  calculated  from  the  PAW  wavefunctions.  (Another  all-electron  DFT 
molecular  dynamics  technique  using  a  mixed-basis  approach  is  applied  in  [84].) 

PAW  is  a  recent  addition  to  the  all-electron  electronic  structure  methods  whose  accuracy  appears  to  be 
similar  to  that  of  the  general  potential  LAPW  approach.  The  implementation  of  the  molecular  dynamics 
formalism  enables  easy  structure  optimization  in  this  method. 


( g)  Illustrative  examples  of  the  electronic  and  optical  properties  of  modem  materials 

As  an  indication  of  the  types  of  information  gleaned  from  all-electron  methods,  we  focus  on  one  recent 
approach,  the  ELAPW  method.  It  has  been  used  to  determine  the  band  structure  and  optical  properties  over 
a  wide  energy  range  for  a  variety  of  crystal  structures  and  chemical  compositions  ranging  from  elementary 
metals  [60]  to  complex  oxides  [85],  layered  dichalcogenides  [86,87]  and  nanoporous  semiconductors  [88]. 
The  k  •  p  formulation  has  also  enabled  calculation  of  the  complex  band  structure  of  the  A1  (100)  surface  [89], 

As  an  illustration  of  the  accuracy  of  the  AFC  ELAPW-fc  • p  method,  we  present  the  dielectric  function  of 
GaAs.  The  dielectric  function  is  a  good  gauge  of  the  quality  of  a  method,  since  not  only  do  the  energies  enter 
the  calculation,  but  also  the  wavefunctions  via  the  matrix  elements  of  the  momentum  operator  — iV.  For  the 
calculation  of  the  dielectric  function  (equation  (Al,3.87))  of  GaAs,  the  conduction  bands  were  rigidly  shifted 
so  that  the  highest  peak  agreed  in  both  experiment  and  theory,  a  shift  of  0.75  eV.  The  imaginary  part  of  the 
dielectric  function  is  shown  in  figure  B3.2.5.  Comparing  the  energy  differences  between  the  three  peaks,  we 
find  that  they  agree  to  within  2  meV.  For  a  wider  comparison,  we  plot  the  results  of  two  more  experiments 
(which  only  have  measured  the  two  peaks  at  lower  photon  energy)  and  several  all-electron  calculations  of  the 
dielectric  function  of  GaAs  in  figure  B3.2.6.  The  FLAPW  results  agree  almost  exactly  with  the  AFC  ELAPW 
values.  The  discrepancies  compared  to  the  experimental  results  found  for  the  other  methods  are  considerably 
larger  than  for  the  general  potential  LAPW  results,  particularly  for  E\ . 

A  recent  study  of  a  class  of  nanoporous  materials,  the  cetineites  [88],  offers  further  illustration  of  the 
possibilities  offered  by  the  modern  band  structure  methods.  The  crystal  is  constructed  of  tubes  of  0.7  nm 
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Figure  B3.2.5.  The  imaginary  part  of  the  dielectric  function  of  GaAs,  according  to  the  AFC  ELAPW-fc  *  p  method  (solid 
curve)  [195]  and  the  experiment  (dashed  curve)  [196].  To  coifect  for  the  band  gap  underestimated  by  the  local  density 
approximation,  the  conduction  bands  have  been  shifted  so  that  the  E2  peaks  agree  in  theory  and  experiment. 


(1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9) 

Figure  B3.2.6.  The  energies  of  the'  Ev  and  E\  peaks  relative  to  the  E2  peak  of  the  imaginary-part  of  the  dielectric  function 
of  GaAs,  calculated  by  self-consistent  DFT  all-electron  methods.  These  energies  do  not  depend  on  the  gap  size.  The 
theoretical  methods  are  noted,  as  are  experimental  results  obtained  by  ellipsometry  (see  chapter  B 1 .26).  The  lower  (upper) 
histogram  gives  the  energy  of  peak  E\  (E\)  relative  to  E2.  LCGO  designates  a  linear-combination-of-Gaussian-orbitals 
method,  OLCAO  an  orthogonalized  linear-combination-of-atomic-orbitals  approach.  Sources:  (1)  [195]’  (2)  [196], 

.  (3)  [197],  (4)  [198],  (5)  [199],  (6)  [200],  (7)  [199],  (8)  [201],  (9)  [202]. 

diameter  arranged  in  a  two-dimensional  hexagonal  structure  with  ‘flattened’  SbSe3  pyramids  arranged  be¬ 
tween  the  tubes  (see  figure.  B3. 2.7).  Cetineites  are  of  potential  technological  interest  because,  singularly 
among  nanoporous  materials,  they  are  semiconductors  rather  than  insulators.  In  figure  B3.2.8,  we  show 
the  comparison  of  the  predicted  density  of  states  to  the  ultraviolet  photoemission  spectrum  (PES,  see  chap¬ 
ter  B  1.1).  The  DOS  can  explain  the  two  main  structures  in  the  PES  at  about  -3  and  -12  eV.  Their  relative 
intensities  agree  with  those  suggested  by  the  DOS  curve.  Three  structures  in  the  DOS  at  -1,  —6  and  -9  eV 
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Cetineite  (Na;Se) 


A6[Sbl20ls][SbX3]2  A= Na,  K;  X=S,  Se 


Figure  B3.2.7.  A  perspective  view  of  the  cetineite  (Na;Se).  The  height  of  the  figure  is  three  lattice  constants  c.  The 
shaded  tube  is  included  only  as  a  guide  to  the  eye.  (From  [88].) 


are  not  resolved  in  the  PES.  This  may  be  due  to  the  selection  rules  of  the  photoemission  process,  not  accounted 
for  in  the  theory,  or  perhaps  due  to  incomplete  angle  integration  experimentally.  The  experimental  results 
confirm,  in  particular,  that  the  number  of  states  is  very  high  close  to  the  valence  band  maximum.  An  orbital 
‘  analysis  shows  that  these  states  are  derived  mainly  from  the  p  states  of  the  O  and  Se  constituents  of  the  crystal, 
with  the  chalcogen  dominating  near  the  top  of  the  valence  band.  Electrons  in  the  Se  p  states  are  thus  most 
easily  excited  into  the  conduction  band.  This,  together  with  their  high  DOS,  makes  the  Se  p  states  located  on 
the  pyramids  the  prime  candidates  for  the  initial  states  of  the  photoconductivity  observed  in  the  cetineites. 

As  another  example  of  properties  extracted  from  all-electron  methods,  figure  B3.2.9  shows  the  results  of 
a  PAW  simulation  of  benzene  molecules  on  a  graphite  surface.  The  study  aimed  to  show  the  extent  to  which 
the  electronic  structure  of  the  molecule  is  modified  by  interaction  with  the  surface,  and  why  the  images  do 
not  reflect  the  molecular  structure.  The  PAW  method  was  used  to  determine  the  structure  of  the  molecule  at 
the  surface,  the  strength  of  the  interaction  between  the  surface  and  the  molecules,  and  to  predict  and  explain 
scanning  tunnelling  microscope  (STM)  images  of  the  molecule  on  the  surface  [90]  (the  STM  is  described  in 
section  B1J9). 

B3.2.3.3  Linear-scaling  electronic  structure  methods 

DFT  calculations  such  as  the  ones  mentioned  in  chapter  A1.3  and  section  B3.2.3.2  become  computationally 
very  expensive  when  the  unit  cell  of  the  interesting  system  becomes  large  and  complex,  with  certain  parts  of 
the  computational  algorithm  typically  scaling  cubically  with  system  size.  A  recent  objective  for  treating  large 
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llTScn  ACpwf°"  °?!  Ph°t0emisS,0n  sPectrum  for  the  cetineite  (Na;Se)  and  the  density  of  states  calculated 
oy  me  Arc  tLAPW-fc  •  p  method  [88]. 


Figure  B2 1.2.9.  A  benzene  molecule  on  a  graphite  surface  [90],  The  geometry  and  the  charge  density  (indicated  by  the 
surfaces  of  constant  density)  have  been  obtained  using  the  PAW  method.  (Figure  by  Professor  P  E  Blochl.) 


sys  ems  is  to  have  the  computational  burden  scale  no  more  than  linearly  with  system  size.  Methods  achieving 
this  are  called  Linear-scaling  or  O(N)  (order  N)  methods,  most  of  which  are  based  on  the  Kohn-Sham  equation 
(see  equation  (B3.2.4)),  aiming  to  calculate  single-electron  wavefunctions,  the  Kohn-Sham  orbitals.  These 
methods  tend  to  be  faster  than  the  conventional  Kohn-Sham  approach  above  a  few  hundred  atoms  [20, 91-93]. 
Another  class  of  methods  is  based  directly  on  the  DFT  of  Hohenberg  and  Kohn  [24],  With  these  techniques  one 
seeks  to  determine  directly  the  density  that  minimizes  the  total  energy;  they  are  often  referred  to  as  orbital-free, 
methods  [94-97],  Such  orbital-free  calculations  do  not  have  the  bottlenecks  present  in  orbital-based  O (N) 
UP  1  calculations,  such  as  the  need  to  localize  orbitals  to  achieve  linear  scaling,  orbital  orthonormalization, 
or  Bnlloum  zone  sampling.  Without  such  bottlenecks,  the  calculations  become  very  inexpensive 
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Equation  (B3.2.3)  lists  the  terms  comprising  the  calculation  of  the  total  energy.  The  term  due  to  the 
external  potential  and  the  Hartree  term  describing  the  Coulomb  repulsion  energy  among  the  electrons  already 
explicitly  depend  on  the  density  instead  of  on  orbitals.  More  difficult  to  evaluate  is  G[p]  =  Ts[p]  +  Exc[p],  a 
functional  which  is  not  known  exactly.  However,  over  the  years  a  number  of  high-quality  exchange-correlation 
functionals  have  been  developed  for  all  kinds  of  systems.  Only  quite  recently  have  more  accurate  kinetic 
energy  density  functionals  (KEDFs)  become  available  [97—99]  that  afford  linear-scaling  computations. 

One  current  limitation  of  orbital -free  DFT  is  that  since  only  the  total  density  is  calculated,  there  is  no  way 
to  identify  contributions  from  electronic  states  of  a  certain  angular  momentum  character  /.  This  identification 
is  exploited  in  non-local  pseudopotentials  so  that  electrons  of  different  l  character  ‘see’  different  potentials, 
considerably  improving  the  quality  of  these  pseudopotentials.  The  orbital-free  methods  thus  are  limited  to 
local  pseudopotentials,  connecting  the  quality  of  their  results  to  the  quality  of  the  available  local  potentials. 
Good  local  pseudopotentials  are  available  for  the  alkali  metals,  the  alkaline  earth  metals  and  aluminium 
[100, 101]  and  methods  exist  for  obtaining  them  for  other  atoms  (see  section  VI.2  of  [97]). 

The  orbital-free  method  has  been  used  for  molecular-dynamics  studies  of  the  formation  of  the  self¬ 
interstitial  defect  in  A1  [102],  pressure-induced  glass-to-crystal  transitions  in  sodium  [103]  and  ion-electron 
correlations  in  liquid  metals  [101].  Calculations  of  densities  for  various  A1  surfaces  have  shown  excellent 
agreement  between  the  charge  densities  as  calculated  by  Kohn-Sham  DFT  and  an  orbital-free  method  using 
a  KEDF  with  a  density-dependent  response  kernel  [99].  The  method  was  used  recently  to  examine  the 
metal-insulator  transition  in  a  two-dimensional  array  of  metal  quantum  dots  [104],  where  the  theory  showed 
that  minute  overlap  of  the  nanoparticle’s  wavefunctions  is  enough  to  transform  the  array  from  an  insulator  to 
a  metal.  As  an  example  of  the  ease  with  which  large  simulations  can  be  performed,  figure  B3.2.10  shows  a 
plot  of  the  charge  density  from  an  orbital-free  calculation  of  a  vacancy  among  255  A1  atoms  [98],  carried  out 
on  a  workstation. 


B3.2.3.4  The  Hartree-Fock  method  in  crystals 

The  HF  method  (discussed  in  section  A  1.3. 1.2)  is  an  alternative  to  DFT  approaches.  It  does  not  include 
electron  correlation  effects*  i.e.  non-classical  electron— electron  interactions  beyond  the  Coulomb  and  exchange 
interactions.  The  neglect  of  these  terms  means  that  the  Coulomb  interaction  is  unscreened,  and  hence  the 
electron  repulsion  energy  is  too  large,  overestimating  ionic  character,  which  leads  to  band  gaps  that  are  too 
large  by  a  factor  of  two  or  more  and  valence  band  widths  that  are  too  wide  by  30-40%  [63].  However,  the 
HF  results  can  be  used  as  a  well  defined  starting  point  for  the  inclusion  of  many-particle  corrections  such  as 
the  GW  approximation  [31,32]  or,  with  considerably  less  computational  effort,  the  results  can  be  improved 
considerably  by  accounting  for  the  Coulomb  hole  and  screening  the  exchange  interaction  using  the  dielectric 
function  [63, 105]. 

Ab  initio  HF  programs  for  crystals  have  been  developed  [106, 1 07]  and  have  been  applied  to  a  wide  variety 
of  bulk  and  surface  systems  [108, 109].  As  an  example,  a  periodic  HF  calculation  using  pseudopotentials 
and  an  LCAO  basis  predicted  binding  energies,  lattice  parameters,  bulk  moduli  and  central-zone  phonon 
frequencies  of  17  III-V  and  IV-IV  semiconductors.  The  authors  find  that  \  .  .[o]n  the  whole,  the  HF  LCAO 
data  appear  no  worse  than  other  ab  initio  results  obtained  with  DF-based  Hamiltonians’  [110].  They  suggest 
that  the  largest  part  of  the  errors  with  respect  to  experiment  is  due  to  correlation  effects  and  to  a  lesser  extent 
due  to  the  imperfections  of  the  pseudopotentials  [1 10].  More  recently,  the  electronic  and  magnetic  properties 
of  transition  metal  oxides  and  halides  such  as  perovskites,  which  had  been  a  problem  earlier,  have  been 
investigated  with  spin-unrestricted  HF  [1 1 1].  In  general,  the  periodic  HF  method  is  best  suited  for  the  study 
of  highly  ionic,  large  band  gap  crystals  because  such  systems  are  the  least  sensitive  to  the  lack  of  electron 
correlation.  "V 


I 'gUre  ®3  '10'  Contour  plot  of  the  electron  densi[y  obtained  ^  an  orbital-free  Hohenberg-Kohn  technique  [98]  The 
figure  shows  a  vacancy  in  bulk  aluminium  in  a  256-site  cell  containing  255  AI  atoms  and  one  empty  site,  the  vacancy 
Dark  areas  represent  low  electron  density  and  light  areas  represent  high  electron  density.  A  Kohn-Sham  calculation  for 

a  cell  of  this  size  would  be  prohibitively  expensive.  Calculations  on  smaller  cell  sizes  using  both  techniques  yielded 
densities  that  were  practically  identical:  H  • 

B3.2.3.5  Quantum  Monte  Carlo 

nv^"  n U6S  p™v*de  accurate  calculations  of  many-electron  systems.  In  variational  QMC  (VMC) 
[112-114],  the  total  energy  of  the  many-electron  system  is  calculated  as  the  expectation  value  of  the  Hamilto- 

!H,amJarameterS  3  mal  waverfunctlon  m  optimized  so  as  to  find  the  lowest-energy  state  (modem  methods 
instead  minimize  the  yanance  of  the  local  energy  [1 15]).  A  Monte  Carlo  (MC)  method  is  used  to  perform 

the  multi-dimensional  integrations  necessary  to  determine  the  expectation  value, 


,  /  JO/pdr 

where  y  is  the  trial  wavefunction  and  |  'P  |2/  / 1  |2dr  is  a  normalized  probability  distribution.  The  integration 
is  pe  orme  y  summing  up  the  local  energy  at  points,  corresponding  to  electron  configurations,  given  by 
fte  probability  distribution.  A  random  walk  algorithm,  such  as  the  Metropolis  algorithm  [116],  is  used  to 
sample  those  regions  of  configuration  space  more  heavily  where  the  probability  density  is  high.  The  standard 
Slater  Jastrow  trial  wavefunction  is  the  product  of  a  Slater  determinant  of  single-electron  orbitals  and  a 
Jastrow  factor,  a  function  which  includes  the  description  of  two-electron  correlation.  As  an  example,  the  trial 
wavefunction  used  for  a  silicon  crystal  contained  32  variational  parameters  whose  optimization  required  the 
calculation  of  the  local  energy  for  10000-20000  statistically  independent  electron  configurations  [117],  In 
conriast  to  the  DMC  technique  described  below,  the  accuracy  of  a  VMC  calculation  depends  on  the  quality 
of  the  many-particle  wavefunction  used  [114].  In  figure  B3.2.11,  we  show  the  determination  of  the  lattice 
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Figure  B3.2.11.  Total  energy  versus  lattice  constant  of  gallium  arsenide  from  a  VMC  calculation  including  256  valence 
electrons  [1 18];  the  curve  is  a  quadratic  fit.  The  error  bars  reflect  the  uncertainties  of  individual  values.  The  experimental 
lattice  constant  is  10.68  au,  the  QMC  result  is  10.69  (±0.1)  au  (Figure  by  Professor  W  Schattke). 


constant  of  GaAs  by  VMC  by  minimization  of  the  total  energy  [118].  This  figure  illustrates  the  roughness  of 
the  potential  energy  surface  due  to  statistical  errors,  which  poses  a  challenge  then  for  the  calculation  of  forces 

with  QMC.  .  . 

In  the  diffusion  QMC  (DMC)  method  [114, 1 19],  the  evolution  of  a  trial  wavefuncfion  (typically  wave- 

functions  of  the  Slater-Jastrow  type,  for  example,  obtained  by  VMC)  proceeds  in  imaginary  time,  r  —  if, 
according  to  the  time-dependent  Schrodinger  equation,  which  then  becomes  a  diffusion  equation.  All  com¬ 
ponents  of  the  wavefuncfion  except  for  the  ground-state  wavefuncfion  are  damped  by  the  time  evolution 
operator  exp(-itfr)  =  exp(-Hr).  The  DMC  was  developed  as  a  simplification  of  the  Green’s  function  MC 
technique  [113]  A  particularly  well  known  use  of  the  Green’s  function  MC  technique  was  the  determination 
by  Ceperley  and  Alder  of  the  energy  of  the  uniform  electron  gas  as  a  function  of  its  density  [120].  This  E(p) 
was  subsequently  parametrized  by  Perdew  and  Zunger  for  the  commonly -used  LDA  exchange-correlation 
potential  [40].  Usually  two  approximations  are  made  to  make  DMC  calculations  tractable:  the  fixed-node 
approximation,  in  which  the  nodes,  the  places  where  the  trial  function  changes  sign,  are  kept  fixed  for  the 
solution  to  enforce  the  fermion  symmetry  of  the  wavefuncfion  and  the  so-called  short-time  approximation, 
whose  effect  can  be  made  very  small  [114],  Excited  states  have  been  calculated  by  replacing  an  orbital  m  the 

Slater  determinant  of  the  trial  wavefuncfion  by  a  conduction-band  orbital  [121].  . 

Recently,  a  method  has  been  proposed  to  overcome  the  problems  associated  with  calculating  forces  in 
both  VMC  and  DMC  [122].  It  has  been  suggested  that  the  use  of  QMC  in  the  near  future  to  tackle  the 
energetics  of  systems  as  challenging  as  liquid  binary  iron  alloys  is  not  unthinkable  [123] . 

B3.2.3.6  Summary  and  comparisons 

As  we  have  outlined,  a  very  wide  variety  of  methods  are  available  to  calculate  the  electronic  structure  of 
solids.  Empirical  TB  methods  (such, as  discussed  in  section  B3.2.2)  are  the  least  expensive,  affording  the 
calculation  of  unit  cells  with  large  numbers  (e.g.  103)  of  atoms,  or  to  provide  cheap  input  to  subsequent 
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methods,  at  the  price  of  quantitative  accuracy.  DFT  methods  (sections  a\*ka  ^ 

10  “  *—  -  ^  as  S:3  ^  “™  *  «*  -MU  ^ 

found  for  DFT.  For  solids,  HF  appears  to  yield  results  JenlnoD^toZT^  “T*  *  **  d liv«% 
but  being  a  genuine  many-particle  theory  it  offers  the  possibilitv  for  r  •  *  glect  of  electron  correlation 
Finally,  the  QMC  techniques  (section  B3.2.3.4)  hold  promise  for  oenn™8^  COrrections- in  contrast  to  DFt’ 
are  still  far  from  able  to  offer  the  same  quanS^^ 

theories  mentioned  before.  With  this  wide  range  of  methods  now  intmd,  ,  enals  and  geometries  as  the 
to  chemisorption  on  solid  surfaces.  traduced,  we  wiU  lo°k  at  their  application 

B3.2.4  Quantum  structural  methods  for  solid  surfaces 

B32A.1  Introduction 

fordeuT  ^  iniH0  qUantUm  dhe^s^n^^^S<^figumionehiteractioii ^Cll^ertmh CU^SS °D dl0Se surfaces 
slabs.  In  betweCenaEeBt3wo  extremes  ^^^1!^  modelf  %  h^ 

the  drastic  approximation  made  by  using  a  finite  cluster  to  describe  re^°gnize  and  attemPt  to  coirect 

electronic  structure  is  inherently  delocalized  or  an  ionic  crystal  with  1/  ®X3mp  e’ a  metallic  conductor  whose 
chemisoiption,  the  binding  of  an  atom  or  a  molernlp  tn  .  n£~ranSe  Coulomb  interactions.  Upon 

in  the  bond  between  the  a^sorbate^ Shari»g  * 
localization  of  the  electrons.  The  attractive  feature  of  the  embedded  f  *6  CrySta  symmetry  wi»  induce 
strengths  of  the  cluster  approach,  namely  it  allows  one  to  describe  the  C  ldea  1S  that  il  Preserves  the 
a  high  degree  of  accuracy  by,  for  example,  quantum  che^S  method,  P^cess  of  chemisorption  to 

account  for  the  presence  of  the  rest  of  the  surface  and  bulk  Surfer  ’  W  6  ^ the  Same  time  attempting  to 
have  been  studied  on  a  variety  of  surfaces,  including  insulators  sefo01^1^1011  a“d  moIecular  adsorption 
these  methods,  we  will  focus  on  those  used  to  examine  adsomfef’  SfemiConductors  and  metals.  To  illustrate 
surfaces.  This  is  not  a  comprehensive  review  oS  appTch"  ^  m°leCUleS  °n  transition  metal 

demonstrate  the  range  of  techniques  and  applicators,  and  some  of  “™Ple!  “ 

B3.2.4.2  The  finite  cluster  model 

clust^mf txansftion'memt afoins^  rMging^rom  vlTalST'  “  *  *"*  adSOrption  °n  3  smaI1’  finite 
calculations  can  be  performed,  typically  the  core  electrons  of  * f0mS  UP  £°  atoms.  Though  all-electron 
core  potential  (ECP,  the  quantum  chemdstry^vers^OT  of  ™et3^  at°ms  are  replaced  by  an  effective 

core-valence  electron  interaction),  while  the  valence  electrons  °fP°teatiaI  that  accounts  approximately  for  the 
a  HF,  Cl,  PT,  or  DFT  formalism.  Typically  a  few  atoms  m  f“Chmetal  at°m  are  treated  explicitly  within 
all)  electrons  explicitly,  while  surrounding  atoms  tend  to  befecribed^011  C°ntain  ^  ValenCe  (or 
one-electron  ECP  representation,  model  pseudopotentials  or  in  the  r  m°/e  CrUdely  Wlth’  for  example,  a 
point  charges.  Generally,  the  structure  of  the  cluster  is  rhr.se’  j,  ^  °f  10niC  crystals’  a  finite  array  of 
of  this  type  of  approach  include  the  early  work  of  Unton  an  drZ/Jf*  fragment  of  the  bulk.  Examples 
electronegative  and  electropositive  atoms  on  a  Ni  rlL  andGoddard  H24],  who  examined  adsorption  of 

In  this  model,  only  the  4s  electrons  on  each  Ni  a°tnm  deSlgned  t0  to1™0  vanous  low-index  faces  of  Ni. 

on  each  Ni  atom  were  treated  explicitly,  while  the  3d  electrons  were 


frequencies  and  binding  eiJ^Ba^r J  STp^blih  P^efeiTed  binding  sites’  °eometries.  vibrational 
it  is  more  accurate  to  teat  mem)  'T***  COmparison  study  showing  that 

is  sufficient  to  describe  the  suirou^  level;  while  it 

cluster  should  be  ‘bond-prepared’,  namely  that  one  sh„S  ?  T  "  [126]  ProPosed  the  idea  that  a 

has  enough  singly-occupied  orbitals  of  the  cornet  S  Udy  ^  e,ectromc  state  of  the  finite  cluster  that 
form  the  necessary  covalent  bonds  between  the  ,t0  interact  with  the  incoming  admolecule  to 

■  surface  reaction,  Panas  «  al  [127  SSSSS?8  T*  6  In  °ne  of  ^  first  stad^  of  a  metal 
level  for  02  on  a Ni,3  cluster,  *  the  Reference  Cl 

used  DFT-LDA  with  a  Gaussian  basis  to  examine  rhe  3  •  bUt  ^  4s  ejfctrons-  Salahub  and  co-workers  [128] 
containing  up  to  16  atoms  mearn  ?ZusZ7Tl°f  Cf^.  “d  HC°°  0n  Nl  <*** 

scheme  improved  dramatically  the  binding  eneroiec  fn  h  a  d  f  .?eS  °f  Nl'  Gradient  corrections  to  the  LDA 
to  experimental  results  for  Ni(l  1 1)  and  nT(  100M 1 29]'  MultelT  ^  Sma11  N*  dUSterS’  When  comPared 
for  example,  in  the  case  of  hydrogen  on  Pd  StaS 

calculated  by  DFT-LDA  for  clusters  contain!™  „n  tn  i-  <=>  ™U1U)  [130],  Diffusion  barriers  were  also 

examples  include  HF  calculations  of  K  adsorbed  on  P.  i  °f  Pd’  Rh’  Sn’  and  Zn  H31].  Other 

PT  (MP2)  calculations  of  acetylene  on  Cu  and  Pd  dusters'  [133  mriTfi  IT*  Mf  ®r_Ples.set  sec^d-order 
calculations  for  CO  on  Cu  clusters  rmi  ^  [  3J’  modified  coupled  pair  functional  (CPF) 

Cu  clusters  [135],  HF,  CASSCF  (complete  active^nar.e  “Culatl0ns  of  hydrogen  adsoiption  on  relaxed 
calculations  for  CO  [136]  and  O  [137]  on  Pt  cluster*  Self'consistent  field>  and  multireference  Cl  and  PT 
tetramers  [138]  and  of  K  and  CO  on  Pds  "  h  W*P*nz«l  ?FT  of  c-CH2N2  on  Pd  and  Cu 

systematically  include  hi°h  levels  of  electron  ro  vantage  of  the  finite-cluster  model  is  that  one  can 

band  structure,  the  presence  ^  “  *  be  baIanCed  against  *e  Iac*  *  a  proper 

.  Next  we  outline  cuirent  strategie^for^meliXating^ome  oVthe^difficulti^^^ t0  modeP*n2  *ow  coverages. 

B3. 2.4.3  Finite -cluster  model  in  contact  with  a  classical  background 
Several  modifications  of  the  finite-cluster  mndAi  mflon f  «. 

and  to  compensate  for  the  lack  of  a  proper  band  strum  °haCC0Unt  for  the  background  Fermi  sea  of  electrons 
mahons  of  the  surface/bulk,  TheJ  approxi- 

crystals  (see,  for  example,  [140]).  Of  these  the  mnHfi  ■  . ,  in  eractlons  and  usually  applied  to  ionic 

sidered  adsoiption  on  metal  surfaces  [141]'  The  so  r  M^vr  \  akatSuji  1S  the  Primary  one  that  has  con- 
cfuSter  pl„s  Z  adsorbate  as  Stet.’to  '  W  nSef "fT  '‘42J  *  “““ 

A  normal  HF  calculation  on  the  small  system  k  nf-rf  a  *  u  ermi  SCa  e*ectrons  °f tbe  bulk  metal. 

■he  cluster  i„  each  “  °'™>°vcd  from 

electron  transfer,  dE/dn,  to  the  work  function  of  the  “  ,  the  total  energy  with  respect  to  the  fractional 
add  us  ter  and  the  bulk  metal  can  be  established  Thus  h  3  ’  ^  eXtent  of  electron  transfer  between  the 

charge  correction  is  also  accounted  for  In  certain  cases  T  a  Sma11  cIuster  optimized  and  an  image 
and  the  ‘surroundings’;  then  ^  “  tninsfcmsd  between  the 

purely  classical  electrostatic  CXample  CI’  can  be  ca^d  out.  This  is  a 

explicit,  manner.  Nakatsuji  has  used  this  to  study  °adsomtionaofSrOUnd  in  an  imPlicit’  rather  than 

one  can  describe  the  polarization  of  the  mem!  10niC  adsorbates  on  metals,  and  finds  that 

[143],  but&u«ia«tC^^^1SZ^  Wm  ^  haVe  WOrked  briefly  With  this  approach 

of  an  ambiguity  of  where  to  pface  ^  tW°-dimensional  cI“ters,  because 

two-dimensional  clusters.  It  fs  also  likely  thaf  he  wavefimr^^f5^*  Suexamples  are  always  smaU  one-  or 
«.oms,  wou,d  not  ad=qM,.,y  ^ 

rder  to  determine  a  cluster  Fermi  level  within  DFT, 
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Figure  B3.2.12.  Schematic  illustration  of  geometries  used  in  the  simulation  of  the  chemisorption  of  a  diatomic  molecule 
on  a  surface  (the  third  dimension  is  suppressed).  The  molecule  is  shown  on  a  surface  simulated  by  (A)  a  semi-infinite 
crystal,  (B)  a  slab  and  an  embedding  region,  (C)  a  slab  with  two-dimensional  periodicity,  (D)  a  slab  in  a  supercell  geometry 
and  (E)  a  cluster. 

originally  by  the  X*  method  (a  simplified  version  of  DFT-LDA;  see  section  Al.3.3.3).  Recent  applications  of 
this  method  have  utilized  more  accurate  forms  of  gradient-corrected  spin-polarized  DFT  to  look  at  adsorption 
of,  for  example,  acetylene  on  Nii4>2o  clusters  [145],  CO  adsorption  on  Ni,  Pd  and  Pt  clusters  of  eight  or  nine 
atoms  [146]  and  NO  adsorption  on  Ru  [147]. 

B3.2.4A  Slab  calculations 

The  other  extreme  of  modelling  chemisorption  is  to  use  a  slab  described  by  DFT  or  HF.  The  slab  is  typically 
taken  to  be  periodic  in  the  directions  parallel  to  the  surface  and  contains  a  few  atomic  layers  in  the  direction 
normal  to  the  surface.  For  the  adatoms  not  to  influence  each  other,  unless  that  is  intended,  the  unit  cell 
needs  to  be  sufficiently  large,  parallel  to  the  surface.  For  computational  reasons,  it  is  advantageous  in  some 
methods,  namely  plane  wave  techniques,  to  have  periodicity  in  three  dimensions.  In  the  supercell  geometry, 
this  periodicity  is  gained  by  considering  slabs  which  are  periodic  in  the  direction  perpendicular  to  the  surface 
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but  separated  from  each  other  by  vacuum  regions.  The  vacuum  region  has  to  be  thick  enough  so  that  there  is 
no  influence  between  the  surfaces  facing  each  other  (the  same  is  true  for  the  slab  thickness).  For  a  schematic 

description  of  several  simulation  model  geometries,  see  figure  B3.2. 12. 

Freeman  and  co-workers  developed  the  FLAPW  method  (see  B3.2.2)  during  the  early  1980s  [70, 1481. 
This  was  a  major  advance,  because  the  conventional  ‘muffin-tin’  potential  was  eliminated  from  their  calcu¬ 
lation  allowing  general-shape  potentials  to  be  evaluated  instead.  Freeman’s  group  first  developed  this  for 
thin  films  and  then  for  bulk  metals.  As  mentioned  before,  the  LAPW  basis,  along  with  the  elimination  of 
any  shape  approximations  in  the  potential,  allows  for  highly  accurate  calculations  on  transition  metal  sur¬ 
faces,  within  the  DFT-LDA  and  the  generalized  gradient  approximation,  GGA  (see  section  A1  3  3  3)  For 
the  ‘stand-alone’  slab  geometry,  figure  B3.2. 12(C),  the  LAPW  basis  functions  decay  exponentially  into  the 
vacuum.  The  numerous  interfacial  systems  examined  by  Freeman’s  group  include,  for  example,  CO  with  K 
or  S  coadsorption  on  Ni(001)  [149],  adsorption  of  sulfur  alone  on  Ni(001)  [150],  Fe  monolayers  on  Ni(l  11) 
[151],  Ag  monolayers  on  MgO(OOl)  [152],  Au-capped  Fe  monolayers  on  MgO(OOl)  [153],  NO  adsorption 
on  Rh,  Pd  and  Pt  [154],  and  Li  on  Ru(001)  [155].  Typical  properties  predicted  are  the  equilibrium  positions, 
magnetic  moments,  charge  densities  and  surface  densities  of  states. 

More  recently,  other  groups  primarily  in  Europe — have  begun  doing  pseudopotential  plane  wave  (often 
gradient-corrected)  DFT  supercell  slab  calculations  (figure  B3.2. 12(D))  for  chemisorption  on  metals.  The 
groups  of  N0rskov  [156],  Scheffler  [157],  Baerends  [158, 159]  and  Hafner  and  Kresse  [160, 161]  have  been 
the  most  active.  Adsorbate-metal  surface  systems  examined  include:  alkalis  and  N2  on  Ru  [156],  NO  on  Pd 
[156],  H2  on  A1  [156],  Cu  [156, 158],  Pd  [157, 158]  and  sulfur-covered  Pd  [157],  CO  oxidation  on  Ru  [157] 

CO  on  Ni,  Pd  and  Pt  [158],  O  on  Pt  [160],  and  H2  on  Rh,  Pd  and  Ag  [161], 

An  interesting  study  by  te  Velde  and  Baerends  [159]  compared  slab-  and  cluster-DFT  results  for  CO 
absorption  on  Cu(100).  They  found  large  oscillations  in  the  chemisorption  binding  energy  of  CO  to  finite 
copper  clusters  as  a  function  of  cluster  size.  This  suggests  that  the  finite-cluster  model  (figure  B3.2. 12(E)) 
is  likely  to  be  inadequate,  at  least  for  modelling  metal  surfaces.  By  contrast,  the  slab  calculations  converge 
quickly  with  the  number  of  Cu  layers  for  the  CO  heat  of  adsorption  and  CO-CO  distances. 

The  supercell  plane  wave  DFT  approach  is  periodic  in  three  dimensions,  which  has  some  disadvantages: 

(l)  thick  vacuum  layers  are  required  so  the  slab  does  not  interact  with  its  images,  (ii)  for  a  tractably  sized  unit 
cell,  only  high  adsorbate  coverages  are  modelled  readily  and  (iii)  one  is  limited  in  accuracy  by  the  form  of  the 
exchange-correlation  functional  chosen.  In  particular,  while  DFT,  especially  using  gradient-corrected  forms 
of  the  exchange-correlation  functional  (GGA),  has  proven  to  be  remarkably  reliable  in  many  instances,  there 
are  a  number  of  examples  for  chemisorption  in  which  the  commonly  used  GGAs  have  been  shown  to  fail 
dramatically  (errors  in  binding  energies  of  ~1  eV  or  greater)  [162, 163].  This  naturally  motivates  the  next  set 
of  approaches,  namely  the  embedded  cluster  strategy. 

B3.2.4.5  Embedded-cluster  schemes:  cluster  in  cluster 

Whitten  and  co-workers  developed  a  metal  cluster  embedding  scheme  appropriate  for  Cl  calculations  during 
the  1980s  [164],  In  essence,  the  method  consists  of:  (i)  solving  for  a  HF  minimum  basis  set  (one  4s 
orbital/atom)  description  of  a  large  cluster  (e.g.,  -30-90  atoms);  (ii)  localizing  the  orbitals  via  exchange 
energy  maximization  with  atomic  basis  functions  on  the  periphery;  (iii)  using  these  localized  orbitals  to  set  up 
effective  Coulomb  and  exchange  operators  for  the  electrons  within  the  cluster  to  be  embedded;  (iv)  improving 
the  basis  set  on  the  atoms  comprising  the  embedded  cluster  and  (v)  performing  a  small  Cl  calculation  (O(103) 
configurations)  within  orbitals  localized  on  the  embedded  cluster.  This  strategy  provides  an  approximate 
way  of  accounting  for  nearby  electrons  outside  the  embedded  cluster  itself.  Whitten  and  co-workers  have 
applied  it  to  a  variety  of  adsorbates  (H,  N,  O,  C — containing  small  molecules)  on,  primarily,  Ni  surfaces 
Duarte  and  Salahub  recently  reported  a  DFT-cluster-in-DFT-cluster  variant  of  Whitten’s  embedding  with  a 
couple  of  twists  on  the  original  approach  (for  example,  fractional  orbital  occupancies  and  charges,  and  an 
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ra,n^en°?aliZinS  the  Ch3rge  0n  the  cluster-  ^  was  implemented  St  ?  §  S  s"miorthogonaJ  basis 

[167  ,  and  by  broadening  each  discrete  energy  level  to  mimic  the  bulk  h^d  ^  method  ln  the  fo™er  case 
the  cluster  in  the  latter  case  [168],  “  buJk  band  structure  within  HF  theory  f0r 

Frsani  [i(59]  has  used  the  density  of  states  from  periodic  HFfseeFn  ?  ?  a\  i  k  , 
the  host  in  which  the  cluster  is  embedded,  where  the  applications  hat  i  3  Slab,caIculations  “>  describe 
as  LiF.  The  onginal  calculation  to  derive  the  external  Coulomb  „  h  ,  b  pnmanly  t0  10nic  crystals  such 

function/density  with  a  proper  band  structure  would  be  preferable^  Seed ^:d“onall>'  “fi™e  wave- 
out  recently  a  series  of  problems  with  such  cluster-in-clus^eremhlnH1  d  '  R°SCh  co'workers  pointed 
of  marked  improvement  of  the  results  over  finite  clusters  of  the  §-  approaches-  These  “elude  the  lack 
partitioning  such  that  charge Conservation rT  **}  ^  ^  °thM  Space 

-tax  [170],  the  inherent  delocalized  nature  ^  ^  . 

B3. 2.4.6  Embedding  of  clusters. in  periodic  background 
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the  near-surface  layers  using  an  embedding  potential  constructed  f*  ^  ®  d’®fnesh-and  co-workers  embed 
all-electron  approach,  using  an  LAPW  basis  [173]  Scheffler  and  V™™  i  u  k  Green’s  function  within  an 

method  GaUSSian  baSis  fDr  the  valence  .electrons  and  pseu^opotenti^Tmi8  approacil 

method  is  somewhat  different  from  In<desfield’s  and  Renechv  •  ,u  *74k  Iormulation  of  the  latter 

die  Green’s  function  and  density  are  known  (typically  the  bulk  metalf  &  ^  Srence  sys,tem  is  chosen  for  which 
in  order  to  get  a  A(embedding  potential)  and  hence  a^fdensitv)  Th’  3  A(Green’s  function)  is  solved  for 
potential  locally  in  a  small  region  around  the  adsorbate  Thesem  Jh^  ^n  W$  T&  t0  Solve  for  the  embedd“g 
calculation  of  the  embedding  density  which  yields  a  'trust  rti?*!^  3  °W  f°r  311  economical  yet  accurate 

dopotentials  to  describe  adsorb  on^Smensbnally  in^T"  “Tl™2  method  using  LDA  with  pseu- 
Williams  er  «/  [176],  The  physical  baZforth^nl I  T  T*  Sbbs  [175k  based  earlier  wo*  by 
off  which  the  Bloch  waves  of  the  perfect  substrate  scatter'5  The  i  f  adsorbate  may  be  considered  a  defect 
of  screening  by  the  electron  gas  of  the  metal  Feibelimn  h  in  eractlon  region  is  short-range  because 
the  chemisorption  of  aaH,  molecule  on  Rh(00n  m;,  e  '1  10  tor 
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surfaces  by  smoothing  out  the  step  function  that  determines  the  occupation  numbers  near  the  Fermi  level. 
Keys  to  the  numerical  success  of  their  method  included:  (i)  symmetric  orthogonalization  of  the  Bloch  basis 
to  produce  a  localized  set  of  functions  that  yielded  a  balanced  distribution  of  charge  in  the  system  and 
(ii)  self-consistent  evaluation  of  the  Fermi  energy  by  fixing  the  charge  on  the  cluster  to  be  neutral.  The  slab 
was  described  with  a  Slater  basis  at  the  DFT-LDA  level,  while  the  embedded  cluster  orbitals  were  expanded  in 
terms  of  Gaussian  functions  at  the  DFT-LDA  level.  While  some  properties  exhibited  non-monotonic  behaviour 
with  increasing  cluster  size,  the  charge  transfer  between  the  metal  surface  and  the  adsorbate  seemed  to  be 
well  described.  They  concluded  that  properties  are  not  well  converged  in  this  method  if  the  cluster  does  not 
contain  shells  of  metal  atoms  that  are  at  least  next-nearest-neighbours  to  the  adsite  metal  atoms. 

Head  and  Silva  used  occupation  numbers  obtained  from  a  periodic  HF  density  matrix  for  the  substrate  to 
define  localized  orbitals  in  the  chemisorption  region,  which  then  defines  a  cluster  subspace  on  which  to  carry 
out  HF  calculations  [181].  Contributions  from  the  surroundings  also  only  come  from  the  bare  slab,  as  in  the 
Green’s  matrix  approach.  Increases  in. computational  power  and  improvements  in  minimization  techniques 
have  made  it  easier  to  obtain  the  electronic  properties  of  adsorbates  by  supercell  slab  techniques,  leading  to 
the  Green’s  function  methods  becoming  less  popular  [182]. 

Cortona  embedded  a  DFT  calculation  in  an  orbital -free  DFT  background  for  ionic  crystals  [183],  which 
necessitates  evaluation  of  kinetic  energy  density  functionals  (KEDFs).  Wesolowski  and  Warshel  [184]  had 
similar  ideas  to  Cortona,  except  they  used  a  frozen  density  background  to  examine  a  solute  in  solution  and 
examined  the  effect  of  varying  the  KEDF.  Stefanovich  and  Truong  also  implemented  Cortona’s  method  with 
a  frozen  density  background  and  applied  it  to,  for  example,  water  adsorption  on  NaCl(OOl)  [185]. 

B3.2.4.7  Embedding  explicit  correlation  methods  in  a  DFT  background 

In  principle,  DFT  calculations  with  an  ideal  exchange-correlation  functional  should  provide  consistently 
accurate  energetics.  The  catch  is,  of  course,  that  the  exact  exchange-correlation  functional  is  not  known. 
While  various  GGAs  have  been  remarkably  successful,  there  are  notable  exceptions  [186, 187],  including 
ones  specific  to  surface  adsorption  mentioned  earlier,  where  the  binding-energy  errors  can  be  more  than  an 
eV  [162, 163].  As  another  example,  Louie  and  Cohen  and  co-workers  found  no  systematic  improvement  over 
the  LDA  when  gradient  corrections  were  included  in  calculations  of  Al,  Nb  and  Pd  bulk  properties,  including 
the  cohesive  energy  [186].  Indeed,  the  design  of  exchange-correlation  functionals  constitutes  an  active  field 
of  research  (see,  for  example,  [188]).  The  lack  of  completely  systematic  means  to  improve  these  functionals 
is  an  unappealing  aspect  of  these  calculations. 

A  first  step  towards  a  systematic  improvement  over  DFT  in  a  local  region  is  the  method  of  Aberenkov  et 
al  [189],  who  calculated  a  correlated  wavefunction  embedded  in  a  DFT  host.  However,  this  is  achieved  using 
an  analytic  embedding  potential  function  fitted  to  DFT  results  on  an  indented  crystal.  One  must  be  cautious 
using  a  bare  indented  crystal  to  represent  the  surroundings,  since  the  density  at  the  surface  of  the  indented 
crystal  will  have  inappropriate  Friedel  oscillations  inside  and  decay  behaviour  at  the  indented  surface  not 
present  in  the  real  crystal. 

We  have  developed  a  different  first-principles  embedding  theory  that  combines  DFT  with  explicit  cor¬ 
relation  methods.  We  sought  to  develop  a  method  for  treating  bulk  or  surface  phases  that  is  more  accurate 
than  current  implementations  of  DFT.  The  idea  is  to  provide  more  accurate  predictions  for  local  energetics,, 
such  as  chemisorption  binding  energies  and  adsorbate  electronic  excitation  energies.  To  achieve  this,  our 
theory  improves  upon  the  DFT  description  of  electron  correlation  in  a  local  region.  This  is  accomplished  by 
an  embedding  theory  that  treats  a  small  region  within  an  accurate  quantum  chemistry  approach  [190, 191], 
which  interacts  with  its  surroundings  via  an  embedding  potential,  uembed(r')*  This  ^embedM  is  derived  from 
a  periodic  DFT  calculation  on  the  total  system.  It  is  expressed  purely  in  terms  of  orbital-free  DFT  (kinetic 
and  potential  energy)  interaction  terms  between  the  embedded  region  and  its  surroundings  a  la  Cortona  and, 
in  particular,  purely  in  terms  of  functionals  of  the  total  density,  ptot,  and  the  density  of  the  embedded  region, 
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A-  We  thus  avoid  construction  of  localized  orbitals  to  describe  the  electrons  in  the  surrounding  environm 
This  is  especially  important  for  metal  surfaces,  where  the  extensive  fc-point  sampling  required  to  get  a  well 
converged  density  makes  localization  impractical  (very  expensive).  This  way  of  expressing  the  embeddi  o 
operator  also  eliminates  problems  that  occur  in  other  forms  of  embedding,  such  as  those  of  matching  condT 
tions  at  the  embedding  boundary,  or  spurious  charge  transfer,  since  the  electrostatic  potential  and  the°densitv 
are  continuous  by  construction.  Its  only  real  disadvantage  is  that  there  is  an  arbitrariness  associated  with  the 
choice  of  Ts.  Development  of  optimal  Ts  functionals  is  an  active  area  of  research  in  our  group  [97-99] 

The  self-consistent  embedding  cycle  proceeds  as  follows.  First,  a  well  converged  density,  Aot  is  calcu 
lated  for  the  extended  metal  surface  in  the  presence  of  an  adsorbate.  This  is  accomplished  within  a  standard" 
pseudopotential  plane  wave  DFT  calculation  (see  chapter  A 1 .3).  Second,  we  partition  the  system  into  the 
region  of  interest  (typically  the  adsorbate  and  neighbouring  metal  atoms  at  or  near  the  surface)  and  its  sur 
roundings  (all  the  other  atoms  in  the  periodic  unit  cell).  The  embedded  region  is  defined  by  the  integral" 
number  of  electrons  and  nuclei  within  that  region  but  not  by  a  particular  physical,  fixed  boundary.  This  allows 
for  the  electron  density  from  the  embedded  region  to  expand  or  contract  variationally  into  the  surrounding 
thus  affording  some  effective  charge  polarization  to  occur  as  needed.  °  ’ 

The  electron  density,  a ,  of  the  embedded  cluster/adsorbate  atoms  is  calculated  using  quantum  chemistry 
methods  (HF,  PT,  multireference  SCF,  or  Cl).  The  initial  step  in  this  iterative  procedure  sets  reml)ed  (r)  to  zero 
since  a  is  needed  in  order  to  calculate  it.  On  subsequent  iterations,  the  third  step  is  to  use  a  and  ptot  to 
calculate  uembedO”),  then  insert  it,  as  a  one-electron  operator  expressed  in  matrix  form  in  the  atomic  orbital 
basis  of  the  adsorbate/cluster,  into  the  quantum  chemistry  calculation  of  step  two.  and  then  a  is  updated  (via 
the  wavefunction).  We  repeatedly  update  uembed(r)  and  then  a  until  full  self-consistency  is  achieved,  with 
fixed  Aot-  In  this  way,  we  variationally  optimize  both  the  quantum  chemistry  wavefunction  and,  implicitly 
the  density  of  the  surroundings,  subject  to  fixed  Aot.  We  tacitly  assume  that  the  DFT-slab  density  for  the  total’ 
system,  Aol,  is  in  fact  a  good  representation  and  does  not  need  to  be  adjusted. 

We  have  shown  that  our  embedding  total  energies  may  be  written  in  terms  of  the  total  energy  obtained 
in  step  one  (the  DFT  total  energy  for  the  entire  system),  plus  a  correction  term,  that  subtracts  out  the  DFT 
energy  in  the  local  region  l  and  adds  back  in  an  ab  initio  total  energy  for  that  same  region, 

'  c1  embed  _  27 DFT  ,  r^ab initio  zrDFTy 
tot  —  ^tot  +  ~  £ t  ). 


Thus,  another  way  to  think  of  the  embedding  is  that  the  ab  initio  treatment  of  region  I  is  correcting  the  DFT 
results  in  the  same  region,  for  the  same  self-consistent  density.  We  expect,  then,  that  such  a  treatment  should 
reduce,  for  example,  the  famous  LDA  overbinding  problem  (LDA  bond  energies  are  generally  significantly 
overestimated).  We  have  indeed  seen  a  smooth  decrease  in  the  LDA  overbinding  as  a  function  of  increasing 
electron  correlation.  We  benchmarked  the  method  against  nearly  exact  calculations  on  a  small  system  and 
then  further  corroborated  it  on  experimentally  well  studied  chemisorption  systems:  GO  on  transition  metal 
surfaces.  Our  binding  energies  are  in  good  agreement  with  nearly  full  configuration  interaction  in  the  former 
and  experimental  adsorbate  binding  energies  in  the  latter.  Very  recently,  we  have  demonstrated  that  excitation 
energies  for  adsorbed  CO  are  dramatically  improved  compared  to  experiment  upon  inclusion  of  the  embedding 
potential  [192],  In  the  future,  we  hope  this  method  will  provide  a  general  means  for  accurate  predictions  of 
the  local  electronic  structure  of  condensed  matter. 


B3.2.5  Outlook 

Computational  solid-state  physics  and  chemistry  are  vibrant  areas  of  research.  The  all-electron  methods  for 
high-accuracy  electronic  structure  calculations  mentioned  in  section  B3.2.3.2  are  in  active  development,  and 
with  PAW,  an  efficient  new  all-electron  method  has  recently  been  introduced.  Ever  more  powerful  computers 
enable  more  detailed  predictions  on  systems  of  increasing  size.  At  the  same  time,  new,  more  complex  materials 
require  methods  that  are  able  to  describe  their  large  unit  ceHs  and  diverse  atomic  make-up.  Here,  the  new 


orbital-free  DFT  method  may  lead  the  way.  More  powerful  techniques  are  also  necessary  for  the  accurate 
treatment  of  surfaces  and  their  interaction  with  atoms  and,  possibly  complex,  molecules.  Combined  with  recent 
progress  in  embedding  theory,  these  developments  make  possible  increasingly  sophisticated  predictions  of 
the  quantum  structural  properties  of  solids  and  solid  surfaces. 
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